Constant coefficients for the equation - Comparison with Central Difference scheme with IC and BCs.

In [145]:
from IPython.core.display import HTML
css_file = './custom.css'
HTML(open(css_file, "r").read())
Out[145]:
$$ \frac {\partial U}{\partial t} = D\frac {\partial^2 U}{\partial x^2} - v\frac {\partial U}{\partial x} + c*U + d $$$$ IC: U(x, 0) = \exp(-x**2) $$

$$ BC: Periodic: U(0, t) = U(1, t) $$ $$ \frac {\partial U(0, t)}{\partial x} = \frac {\partial U(1, t)}{\partial t} $$

In [130]:
from scipy.linalg import solve
# Define the equation parameters
from scipy.special import erfc
import pandas as pd
import numpy as np
import scipy.integrate
from numpy import exp
import matplotlib.pyplot as plt
In [ ]:
 
In [ ]:
 
In [131]:
def testAB_PS(a, b, c, d, nx, nt):

    L = np.pi
    T = 4.0


    x = np.linspace(-L, L, nx+1)
    t = np.linspace(0, T, nt+1)
    
    u = np.exp(-x**2)#np.sin(np.pi*xin)
    ui = u;
    U_CDI = np.zeros((nx+1, nt+1))
    error_2norm_CDI_PSTD = np.zeros(nt+1)

    # Discretization of space and time
    dx = 2*L/(nx+1)
    # dt = 0.1*dx**2/D#T/(nt-1)
    # dx = L/(nx+1)
    dt = T/(nt+1)


    #plt.figure()
    #plt.plot(u, label='Initial')
    def numerical_scheme_CDI(u_CDI, a, b, c, d, dx, dt):
        alpha = a*dt/dx
        beta = b*dt/dx**2
        nx = len(u_CDI)-1;
        A = np.zeros((nx+1, nx+1))
        bi = d*dt+u_CDI;
        for j in range(1, len(u_CDI)-1):

            A[j, j-1] = -(beta-alpha/2)
            A[j, j] = 1+2*beta-c*dt
            A[j, j+1] = -(beta+alpha/2)

        # Applying first order Periodic Boundary Condition
        A[0, 0]=(1+2*beta);
        A[0, 1]= (-alpha/2-beta)
        A[0, nx]= (alpha/2-beta)    


        A[nx, 0]= (-alpha/2-beta) ;
        A[nx, nx-1]= (alpha/2 - beta)
        A[nx, nx]= 1+2*beta    


        # Applying second order periodic boundary conditions
        # bt = d*dt+u_CDI;
        # kk = np.zeros(2)
        # bi = [*kk, *bt]

        # A[0, 0]=0;
        # A[0, 1]= 1
        # A[0, N]= -1
        # A[0, N+1]=0

        # A[1, 0] = -1
        # A[1,1] = 0
        # A[1, 2] = 1
        # A[1, N-1] = 1
        # A[1, N] = 0 
        # A[1, N+1] = -1


        # Boundary condition
        # A[0, 0] = -3
        # A[0, 1] = 4
        # A[0, 2] = -1
        # bi[0] = 0

        # A[N, N] = 3
        # A[N, N-1] = -4
        # A[N, N-2] = 1
        # bi[N] = 0


        # A[0,0] = 1.0
        # # A[N,N] = 1.0

        # bi[0] = l
        # bi[N] = r


        u_CDI = solve(A, bi)

        return u_CDI



    # Apply the numerical scheme
    # u = np.exp(-((x+0.5)**2)/(0.00125))
    U_CDI = np.zeros((nx+1, nt+1))
    error_2norm_QI_PSTD = np.zeros(nt)

    u_CDI=u;
    plt.plot(x,u)
    for n in range(nt):
        u_CDI = numerical_scheme_CDI(u_CDI, a, b, c, d, dx, dt)
        U_CDI[:, n] = u_CDI

    #plt.plot(u_CDI, label='Final')
    #plt.show()
    #plt.close()

    plt.rcParams['font.size'] = 18
    fig = plt.figure() 
    plt.imshow(U_CDI, cmap='viridis', extent=[0, T, 0, 2*L], aspect='auto')
    plt.colorbar()
    title = "CDI - nx = "+str(nx+1)
    plt.title(title)
    plt.xlabel('Time')
    plt.ylabel('Position')
    plt.show()


    k = 2*np.pi*np.fft.rfftfreq(nx+1, 2*L/(nx))
    k2=k**2;

    # Defining variables
    U_ps = np.zeros((nx+1, nt+1))
    U_psk = np.zeros((nx+1, nt+1))
    #plt.figure()
    #plt.plot(xin, u, color = 'b', label = "initial")
    # left = (0.025/(np.sqrt(0.000625 + 0.02*t)))*(np.exp((-(0.5 - t)**2)/(0.00125 + 0.04*t)))
    # right = (0.025/(np.sqrt(0.000625 + 0.02*t)))*(np.exp((-(1.5 - t)**2)/(0.00125 + 0.04*t)))

    # Solving over time
    for i in range(nt+1): # 
        
        uk = np.fft.rfft(u)

        uk[:] = (uk[:])/(1-(-b*k**2 + a*1j*k+c)*dt) 
        # print(np.shape(uk))    
        u = np.fft.irfft(uk) #+ d*dt
        U_ps[:, i] = u
        error =  U_CDI[:, i] - u
        error_2norm_CDI_PSTD[i] =np.linalg.norm(error)/np.sqrt(nx)


    fig = plt.figure()  
    plt.imshow((U_ps), cmap='viridis',  extent=[0, T, 0, 2*L], aspect='auto')
    plt.colorbar()
    plt.title("CD Scheme with PSTD")
    plt.xlabel('Time')
    plt.ylabel('Position')
    title = "PSTDI - nx = "+str(nx+1)
    plt.title(title)
    plt.grid(True)
    plt.show()
    plt.close()
    
    U_psk = np.zeros((nx+1, nt+1))
    error_2norm_PSTDI_PSTDK = np.zeros(nt+1)
    
    error_2norm_CDI_PSTDK = np.zeros(nt+1)
    # Spatial grid
    m=nx+1                            # Number of grid points in space
    L = np.pi                   # Width of spatial domain
    x = np.arange(-m/2,m/2)*(2*L/m)   # Grid points
    dx = x[1]-x[0]                  # Grid spacing

    # Temporal grid
    tmax=T     # Final time
    N = nt+1       # number grid points in time
    dt = tmax/N   # interval between output times

    xi = np.fft.fftfreq(m)*m*np.pi/L  # Wavenumber "grid"
    xi = 2*np.pi*np.fft.fftfreq(nx+1, 2*L/(nx))
    
    # k = 2*np.pi*np.fft.rfftfreq(nx+1, 2*L/(nx))
    
    
    # (this is the order in which numpy's FFT gives the frequencies)

    # Initial data
    u = np.exp(-x**2)# np.sin(2*x)**2 * (x<-L/4)
    uhat0 = np.fft.fft(u)
    plt.plot(x,u)
    epsilon=b  # Diffusion coefficient
    # a = 1.0       # Advection coefficient

    print(np.shape(xi))
    print(np.shape(u))
    # Now we solve the problem
    for i in range(0,N):
        t = i*dt
        uhat = np.exp((1.j*xi*a - epsilon*xi**2)*t) * uhat0
        u = np.real(np.fft.ifft(uhat)) #+ d*t
        
        U_psk[:, i] = u
        error =  U_CDI[:, i] - u
        error_2norm_CDI_PSTDK[i] =np.linalg.norm(error)/np.sqrt(nx)
        error =  U_ps[:, i] - u
        error_2norm_PSTDI_PSTDK[i] =np.linalg.norm(error)/np.sqrt(nx)

    
    fig = plt.figure()  
    plt.imshow((U_psk), cmap='viridis',  extent=[0, T, 0, 2*L], aspect='auto')
    plt.colorbar()
    plt.title("CD Scheme with PSTDK")
    plt.xlabel('Time')
    plt.ylabel('Position')
    title = "PSTDK - nx = "+str(nx+1)
    plt.title(title)
    plt.grid(True)
    plt.show()
    plt.close()
    
    plt.figure();
    # plt.plot(x, ui, '-', color='g', label = 'Initial'); 
    plt.plot(x, U_CDI[:, int(nt/2)], '-*', color='b', label = 'CDI Method'); 
    plt.plot(x, U_ps[:, int(nt/2)], '-+', color = 'r', label = 'Spectral Method');
    plt.plot(x, U_psk[:, int(nt/2)], '-.', color = 'g', label = 'Ketchinson Spectral Method');
    
    res =  [U_CDI[:, nt-1], U_ps[:, nt-1], U_psk[:, nt-1]]
    
    plt.legend(loc="upper right")
    title = "Comparison - nx = "+str(nx+1)
    plt.title(title)
    plt.show()
    plt.close()
    plt.figure()
    plt.plot(error_2norm_CDI_PSTD[0:nt], '.', label = 'PSTDI_CDI')
    
    plt.plot(error_2norm_CDI_PSTDK[0:nt], '.', label = 'PSTDK_CDI')
    
    plt.plot(error_2norm_PSTDI_PSTDK[0:nt], '.', label = 'PSTDI_PSTDK')
    plt.legend()
    plt.show()
    plt.close()
    return res

Pe = 100 without any source terms

Pe = 100

In [132]:
a = 1.0
b = 0.01
c = 0.0; 
d = 0.0

testAB_PS(a, b, c, d, 63, 1000)
testAB_PS(a, b, c, d, 127, 1000)
res = testAB_PS(a, b, c, d, 255, 1000)

plt.figure();
plt.plot(res[0], '-*', color='b', label = 'CDI Method'); 
plt.plot(res[1], '-+', color = 'r', label = 'Spectral Method');
plt.plot(res[2], '-.', color = 'g', label = 'Ketchinson Spectral Method');
plt.legend(loc="upper right")
title = "Comparison - nx = "+str(255+1)
plt.title(title)
plt.show()
plt.close()
(64,)
(64,)
(128,)
(128,)
(256,)
(256,)
In [ ]:
 

Pe = 1

In [133]:
a = 1.0
b = 1.0
c = 0.0; 
d = 0.0

testAB_PS(a, b, c, d, 63, 1000)
testAB_PS(a, b, c, d, 127, 1000)
testAB_PS(a, b, c, d, 255, 1000)
(64,)
(64,)
(128,)
(128,)
(256,)
(256,)
Out[133]:
[array([0.28631191, 0.28616154, 0.28600805, 0.28585155, 0.28569211,
        0.28552985, 0.28536484, 0.28519721, 0.28502704, 0.28485444,
        0.28467951, 0.28450236, 0.28432309, 0.28414182, 0.28395865,
        0.28377369, 0.28358706, 0.28339886, 0.28320921, 0.28301822,
        0.28282601, 0.2826327 , 0.2824384 , 0.28224322, 0.28204729,
        0.28185073, 0.28165364, 0.28145616, 0.2812584 , 0.28106047,
        0.2808625 , 0.28066461, 0.28046692, 0.28026954, 0.2800726 ,
        0.27987621, 0.27968049, 0.27948556, 0.27929154, 0.27909854,
        0.27890668, 0.27871607, 0.27852684, 0.27833909, 0.27815293,
        0.27796849, 0.27778587, 0.27760517, 0.27742652, 0.27725001,
        0.27707576, 0.27690386, 0.27673443, 0.27656756, 0.27640335,
        0.27624191, 0.27608333, 0.2759277 , 0.27577512, 0.27562569,
        0.27547948, 0.2753366 , 0.27519712, 0.27506113, 0.27492872,
        0.27479995, 0.27467491, 0.27455368, 0.27443633, 0.27432292,
        0.27421353, 0.27410822, 0.27400706, 0.2739101 , 0.27381741,
        0.27372904, 0.27364505, 0.27356547, 0.27349038, 0.27341979,
        0.27335377, 0.27329235, 0.27323557, 0.27318346, 0.27313605,
        0.27309337, 0.27305545, 0.27302231, 0.27299396, 0.27297044,
        0.27295174, 0.27293789, 0.27292889, 0.27292475, 0.27292546,
        0.27293103, 0.27294146, 0.27295673, 0.27297684, 0.27300178,
        0.27303154, 0.27306608, 0.2731054 , 0.27314947, 0.27319827,
        0.27325175, 0.2733099 , 0.27337268, 0.27344004, 0.27351195,
        0.27358837, 0.27366924, 0.27375453, 0.27384417, 0.27393812,
        0.27403632, 0.2741387 , 0.27424522, 0.27435579, 0.27447037,
        0.27458887, 0.27471122, 0.27483736, 0.27496721, 0.27510068,
        0.2752377 , 0.27537819, 0.27552206, 0.27566922, 0.27581958,
        0.27597306, 0.27612957, 0.276289  , 0.27645126, 0.27661626,
        0.27678389, 0.27695406, 0.27712666, 0.27730159, 0.27747873,
        0.277658  , 0.27783927, 0.27802244, 0.2782074 , 0.27839403,
        0.27858223, 0.27877188, 0.27896286, 0.27915507, 0.27934838,
        0.27954268, 0.27973785, 0.27993378, 0.28013035, 0.28032743,
        0.28052491, 0.28072268, 0.2809206 , 0.28111857, 0.28131646,
        0.28151415, 0.28171153, 0.28190847, 0.28210486, 0.28230058,
        0.28249551, 0.28268954, 0.28288254, 0.2830744 , 0.283265  ,
        0.28345424, 0.28364199, 0.28382815, 0.28401259, 0.28419522,
        0.28437591, 0.28455457, 0.28473108, 0.28490533, 0.28507723,
        0.28524667, 0.28541354, 0.28557775, 0.2857392 , 0.28589778,
        0.28605341, 0.28620599, 0.28635543, 0.28650163, 0.28664452,
        0.286784  , 0.28691999, 0.28705241, 0.28718118, 0.28730622,
        0.28742745, 0.28754481, 0.28765822, 0.28776761, 0.28787292,
        0.28797409, 0.28807104, 0.28816374, 0.28825211, 0.28833611,
        0.28841568, 0.28849078, 0.28856137, 0.28862739, 0.28868881,
        0.2887456 , 0.28879771, 0.28884512, 0.2888878 , 0.28892572,
        0.28895887, 0.28898721, 0.28901073, 0.28902943, 0.28904328,
        0.28905228, 0.28905643, 0.28905572, 0.28905015, 0.28903972,
        0.28902445, 0.28900433, 0.28897939, 0.28894964, 0.28891509,
        0.28887577, 0.2888317 , 0.2887829 , 0.28872942, 0.28867127,
        0.28860849, 0.28854112, 0.28846921, 0.28839279, 0.28831192,
        0.28822663, 0.28813698, 0.28804303, 0.28794483, 0.28784245,
        0.28773593, 0.28762535, 0.28751078, 0.28739228, 0.28726992,
        0.28714378, 0.28701393, 0.28688045, 0.28674343, 0.28660294,
        0.28645907]),
 array([0.28657219, 0.28641911, 0.28626277, 0.28610325, 0.28594065,
        0.28577507, 0.2856066 , 0.28543536, 0.28526144, 0.28508494,
        0.28490598, 0.28472466, 0.28454109, 0.28435539, 0.28416765,
        0.28397801, 0.28378656, 0.28359343, 0.28339873, 0.28320258,
        0.2830051 , 0.2828064 , 0.28260662, 0.28240585, 0.28220424,
        0.28200189, 0.28179894, 0.2815955 , 0.28139169, 0.28118764,
        0.28098348, 0.28077931, 0.28057528, 0.28037149, 0.28016808,
        0.27996517, 0.27976287, 0.27956131, 0.27936061, 0.2791609 ,
        0.27896228, 0.27876489, 0.27856884, 0.27837425, 0.27818123,
        0.27798991, 0.27780039, 0.27761279, 0.27742723, 0.27724382,
        0.27706266, 0.27688387, 0.27670756, 0.27653382, 0.27636277,
        0.2761945 , 0.27602912, 0.27586674, 0.27570744, 0.27555132,
        0.27539848, 0.275249  , 0.27510299, 0.27496052, 0.27482168,
        0.27468656, 0.27455524, 0.27442779, 0.2743043 , 0.27418483,
        0.27406946, 0.27395827, 0.27385131, 0.27374864, 0.27365035,
        0.27355647, 0.27346707, 0.2733822 , 0.27330192, 0.27322626,
        0.27315529, 0.27308903, 0.27302753, 0.27297083, 0.27291896,
        0.27287195, 0.27282983, 0.27279263, 0.27276037, 0.27273306,
        0.27271073, 0.27269338, 0.27268104, 0.2726737 , 0.27267136,
        0.27267404, 0.27268173, 0.27269442, 0.27271211, 0.27273479,
        0.27276244, 0.27279505, 0.27283259, 0.27287505, 0.27292239,
        0.2729746 , 0.27303163, 0.27309346, 0.27316004, 0.27323134,
        0.27330732, 0.27338792, 0.27347311, 0.27356282, 0.273657  ,
        0.2737556 , 0.27385857, 0.27396582, 0.27407731, 0.27419296,
        0.27431271, 0.27443648, 0.2745642 , 0.27469578, 0.27483116,
        0.27497025, 0.27511297, 0.27525923, 0.27540894, 0.27556201,
        0.27571835, 0.27587787, 0.27604046, 0.27620604, 0.2763745 ,
        0.27654575, 0.27671966, 0.27689616, 0.27707511, 0.27725643,
        0.27744   , 0.2776257 , 0.27781343, 0.27800308, 0.27819452,
        0.27838765, 0.27858235, 0.27877849, 0.27897597, 0.27917467,
        0.27937445, 0.27957522, 0.27977683, 0.27997917, 0.28018213,
        0.28038557, 0.28058937, 0.28079342, 0.28099759, 0.28120175,
        0.28140578, 0.28160957, 0.28181298, 0.2820159 , 0.2822182 ,
        0.28241976, 0.28262045, 0.28282017, 0.28301879, 0.28321618,
        0.28341223, 0.28360683, 0.28379984, 0.28399117, 0.28418069,
        0.28436829, 0.28455385, 0.28473726, 0.28491842, 0.28509722,
        0.28527354, 0.28544728, 0.28561833, 0.2857866 , 0.28595198,
        0.28611437, 0.28627367, 0.28642979, 0.28658264, 0.28673211,
        0.28687813, 0.2870206 , 0.28715944, 0.28729457, 0.28742589,
        0.28755334, 0.28767684, 0.28779631, 0.28791168, 0.28802288,
        0.28812984, 0.2882325 , 0.28833081, 0.28842469, 0.28851409,
        0.28859896, 0.28867925, 0.2887549 , 0.28882588, 0.28889214,
        0.28895364, 0.28901034, 0.28906221, 0.28910922, 0.28915134,
        0.28918854, 0.28922081, 0.28924812, 0.28927045, 0.2892878 ,
        0.28930014, 0.28930749, 0.28930982, 0.28930714, 0.28929945,
        0.28928676, 0.28926907, 0.28924639, 0.28921874, 0.28918613,
        0.28914859, 0.28910613, 0.28905879, 0.28900658, 0.28894954,
        0.28888772, 0.28882113, 0.28874983, 0.28867385, 0.28859324,
        0.28850806, 0.28841835, 0.28832416, 0.28822555, 0.28812259,
        0.28801533, 0.28790384, 0.28778819, 0.28766844, 0.28754466,
        0.28741694, 0.28728535, 0.28714997, 0.28701088, 0.28686816,
        0.2867219 ]),
 array([0.28771039, 0.28755647, 0.28739927, 0.28723886, 0.28707535,
        0.28690884, 0.28673943, 0.28656723, 0.28639232, 0.28621483,
        0.28603485, 0.2858525 , 0.28566788, 0.28548111, 0.2852923 ,
        0.28510156, 0.284909  , 0.28471476, 0.28451893, 0.28432164,
        0.284123  , 0.28392315, 0.28372219, 0.28352025, 0.28331744,
        0.28311391, 0.28290975, 0.2827051 , 0.28250009, 0.28229483,
        0.28208944, 0.28188406, 0.28167881, 0.2814738 , 0.28126916,
        0.28106502, 0.28086151, 0.28065873, 0.28045681, 0.28025588,
        0.28005606, 0.27985746, 0.27966021, 0.27946443, 0.27927022,
        0.27907772, 0.27888703, 0.27869828, 0.27851157, 0.27832701,
        0.27814472, 0.27796482, 0.27778739, 0.27761256, 0.27744043,
        0.2772711 , 0.27710468, 0.27694126, 0.27678094, 0.27662382,
        0.27647   , 0.27631956, 0.2761726 , 0.27602921, 0.27588947,
        0.27575346, 0.27562127, 0.27549298, 0.27536867, 0.2752484 ,
        0.27513226, 0.27502031, 0.27491262, 0.27480926, 0.27471028,
        0.27461575, 0.27452572, 0.27444025, 0.27435939, 0.27428319,
        0.27421169, 0.27414493, 0.27408297, 0.27402583, 0.27397355,
        0.27392616, 0.27388368, 0.27384616, 0.2738136 , 0.27378602,
        0.27376345, 0.2737459 , 0.27373337, 0.27372588, 0.27372343,
        0.27372602, 0.27373365, 0.27374631, 0.273764  , 0.27378671,
        0.27381442, 0.27384712, 0.27388478, 0.27392739, 0.27397491,
        0.27402733, 0.2740846 , 0.2741467 , 0.27421358, 0.27428521,
        0.27436154, 0.27444253, 0.27452812, 0.27461828, 0.27471293,
        0.27481203, 0.27491551, 0.27502332, 0.27513539, 0.27525164,
        0.27537202, 0.27549644, 0.27562484, 0.27575713, 0.27589324,
        0.27603309, 0.27617658, 0.27632364, 0.27647417, 0.27662808,
        0.27678529, 0.27694569, 0.2771092 , 0.2772757 , 0.27744511,
        0.27761732, 0.27779222, 0.27796971, 0.27814968, 0.27833203,
        0.27851665, 0.27870342, 0.27889223, 0.27908297, 0.27927552,
        0.27946977, 0.27966559, 0.27986288, 0.28006151, 0.28026137,
        0.28046233, 0.28066427, 0.28086707, 0.28107061, 0.28127476,
        0.28147941, 0.28168442, 0.28188968, 0.28209507, 0.28230045,
        0.28250571, 0.28271071, 0.28291535, 0.28311949, 0.28332301,
        0.28352579, 0.2837277 , 0.28392863, 0.28412846, 0.28432706,
        0.28452431, 0.2847201 , 0.2849143 , 0.2851068 , 0.28529749,
        0.28548625, 0.28567296, 0.28585752, 0.28603981, 0.28621972,
        0.28639715, 0.28657198, 0.28674411, 0.28691344, 0.28707987,
        0.28724329, 0.28740361, 0.28756073, 0.28771456, 0.287865  ,
        0.28801196, 0.28815536, 0.28829511, 0.28843111, 0.2885633 ,
        0.2886916 , 0.28881591, 0.28893618, 0.28905233, 0.28916428,
        0.28927197, 0.28937534, 0.28947432, 0.28956885, 0.28965888,
        0.28974435, 0.28982522, 0.28990142, 0.28997292, 0.29003967,
        0.29010164, 0.29015878, 0.29021107, 0.29025846, 0.29030093,
        0.29033846, 0.29037102, 0.2903986 , 0.29042117, 0.29043872,
        0.29045125, 0.29045874, 0.29046119, 0.2904586 , 0.29045097,
        0.29043831, 0.29042062, 0.29039791, 0.2903702 , 0.2903375 ,
        0.29029983, 0.29025723, 0.2902097 , 0.29015728, 0.29010001,
        0.29003791, 0.28997103, 0.2898994 , 0.28982306, 0.28974207,
        0.28965648, 0.28956632, 0.28947167, 0.28937256, 0.28926908,
        0.28916127, 0.2890492 , 0.28893294, 0.28881256, 0.28868814,
        0.28855974, 0.28842744, 0.28829133, 0.28815148, 0.28800799,
        0.28786093])]

Pe = 0.1

In [134]:
a = 0.1
b = 1.0
c = 0.0; 
d = 0.0

testAB_PS(a, b, c, d, 63, 1000)
testAB_PS(a, b, c, d, 127, 1000)
testAB_PS(a, b, c, d, 255, 1000)
(64,)
(64,)
(128,)
(128,)
(256,)
(256,)
Out[134]:
[array([0.27353551, 0.27361733, 0.27370358, 0.27379423, 0.27388921,
        0.27398847, 0.27409195, 0.27419958, 0.2743113 , 0.27442705,
        0.27454675, 0.27467033, 0.27479771, 0.27492883, 0.2750636 ,
        0.27520194, 0.27534377, 0.275489  , 0.27563754, 0.27578931,
        0.2759442 , 0.27610214, 0.27626303, 0.27642676, 0.27659324,
        0.27676237, 0.27693405, 0.27710817, 0.27728462, 0.27746332,
        0.27764413, 0.27782696, 0.2780117 , 0.27819823, 0.27838644,
        0.27857623, 0.27876746, 0.27896004, 0.27915384, 0.27934874,
        0.27954464, 0.2797414 , 0.27993892, 0.28013707, 0.28033574,
        0.2805348 , 0.28073413, 0.28093362, 0.28113314, 0.28133258,
        0.28153181, 0.28173072, 0.28192918, 0.28212707, 0.28232428,
        0.28252069, 0.28271617, 0.28291062, 0.28310391, 0.28329592,
        0.28348655, 0.28367568, 0.28386318, 0.28404896, 0.28423289,
        0.28441487, 0.28459479, 0.28477254, 0.28494801, 0.28512109,
        0.28529169, 0.2854597 , 0.28562501, 0.28578754, 0.28594717,
        0.28610382, 0.28625738, 0.28640778, 0.28655491, 0.28669869,
        0.28683903, 0.28697585, 0.28710907, 0.2872386 , 0.28736436,
        0.28748629, 0.2876043 , 0.28771833, 0.2878283 , 0.28793416,
        0.28803583, 0.28813327, 0.28822639, 0.28831516, 0.28839952,
        0.28847942, 0.2885548 , 0.28862563, 0.28869185, 0.28875344,
        0.28881036, 0.28886256, 0.28891002, 0.28895271, 0.2889906 ,
        0.28902368, 0.28905191, 0.2890753 , 0.28909381, 0.28910743,
        0.28911617, 0.28912002, 0.28911897, 0.28911302, 0.28910218,
        0.28908645, 0.28906585, 0.28904038, 0.28901006, 0.28897492,
        0.28893496, 0.28889022, 0.28884072, 0.28878649, 0.28872756,
        0.28866397, 0.28859577, 0.28852298, 0.28844565, 0.28836383,
        0.28827757, 0.28818692, 0.28809194, 0.28799268, 0.2878892 ,
        0.28778157, 0.28766984, 0.28755409, 0.28743439, 0.28731081,
        0.28718342, 0.2870523 , 0.28691752, 0.28677918, 0.28663735,
        0.28649212, 0.28634358, 0.28619181, 0.28603691, 0.28587896,
        0.28571808, 0.28555434, 0.28538786, 0.28521873, 0.28504705,
        0.28487293, 0.28469647, 0.28451777, 0.28433696, 0.28415412,
        0.28396938, 0.28378285, 0.28359464, 0.28340485, 0.28321361,
        0.28302104, 0.28282724, 0.28263233, 0.28243644, 0.28223967,
        0.28204215, 0.281844  , 0.28164533, 0.28144627, 0.28124694,
        0.28104745, 0.28084793, 0.28064849, 0.28044926, 0.28025035,
        0.28005189, 0.279854  , 0.27965679, 0.27946038, 0.2792649 ,
        0.27907046, 0.27887717, 0.27868515, 0.27849453, 0.2783054 ,
        0.2781179 , 0.27793213, 0.27774819, 0.27756621, 0.2773863 ,
        0.27720855, 0.27703309, 0.27686   , 0.27668941, 0.2765214 ,
        0.27635609, 0.27619357, 0.27603394, 0.27587729, 0.27572373,
        0.27557334, 0.27542621, 0.27528243, 0.27514209, 0.27500527,
        0.27487206, 0.27474254, 0.27461678, 0.27449485, 0.27437684,
        0.27426282, 0.27415284, 0.27404699, 0.27394532, 0.27384789,
        0.27375476, 0.27366599, 0.27358164, 0.27350175, 0.27342636,
        0.27335554, 0.27328931, 0.27322772, 0.27317081, 0.27311861,
        0.27307115, 0.27302846, 0.27299057, 0.2729575 , 0.27292926,
        0.27290588, 0.27288737, 0.27287374, 0.272865  , 0.27286116,
        0.27286221, 0.27286816, 0.272879  , 0.27289472, 0.27291533,
        0.27294079, 0.27297111, 0.27300626, 0.27304621, 0.27309095,
        0.27314045, 0.27319468, 0.27325361, 0.27331719, 0.2733854 ,
        0.27345819]),
 array([0.27329704, 0.27338114, 0.27346981, 0.27356302, 0.27366071,
        0.27376281, 0.27386926, 0.27398   , 0.27409496, 0.27421408,
        0.27433728, 0.27446449, 0.27459563, 0.27473062, 0.27486938,
        0.27501183, 0.27515788, 0.27530744, 0.27546043, 0.27561675,
        0.2757763 , 0.275939  , 0.27610474, 0.27627342, 0.27644494,
        0.27661921, 0.2767961 , 0.27697552, 0.27715736, 0.27734151,
        0.27752786, 0.27771629, 0.2779067 , 0.27809896, 0.27829296,
        0.27848859, 0.27868573, 0.27888425, 0.27908405, 0.27928499,
        0.27948696, 0.27968984, 0.27989349, 0.28009781, 0.28030267,
        0.28050794, 0.28071351, 0.28091924, 0.28112501, 0.2813307 ,
        0.28153619, 0.28174135, 0.28194605, 0.28215019, 0.28235362,
        0.28255623, 0.2827579 , 0.2829585 , 0.28315792, 0.28335603,
        0.28355272, 0.28374787, 0.28394135, 0.28413306, 0.28432287,
        0.28451068, 0.28469636, 0.28487982, 0.28506093, 0.28523959,
        0.28541569, 0.28558912, 0.28575979, 0.28592758, 0.2860924 ,
        0.28625414, 0.28641271, 0.28656802, 0.28671997, 0.28686847,
        0.28701342, 0.28715475, 0.28729237, 0.28742619, 0.28755613,
        0.28768212, 0.28780408, 0.28792193, 0.28803561, 0.28814504,
        0.28825017, 0.28835092, 0.28844724, 0.28853906, 0.28862634,
        0.28870902, 0.28878705, 0.28886039, 0.28892898, 0.28899279,
        0.28905178, 0.28910592, 0.28915517, 0.2891995 , 0.28923888,
        0.2892733 , 0.28930272, 0.28932714, 0.28934654, 0.2893609 ,
        0.28937023, 0.2893745 , 0.28937373, 0.2893679 , 0.28935703,
        0.28934112, 0.28932018, 0.28929422, 0.28926326, 0.28922731,
        0.28918641, 0.28914057, 0.28908981, 0.28903418, 0.28897371,
        0.28890842, 0.28883837, 0.28876359, 0.28868413, 0.28860003,
        0.28851135, 0.28841814, 0.28832045, 0.28821835, 0.28811189,
        0.28800115, 0.28788618, 0.28776706, 0.28764386, 0.28751665,
        0.2873855 , 0.28725051, 0.28711175, 0.28696929, 0.28682324,
        0.28667367, 0.28652068, 0.28636436, 0.28620481, 0.28604211,
        0.28587636, 0.28570768, 0.28553615, 0.28536189, 0.28518499,
        0.28500557, 0.28482373, 0.28463957, 0.28445323, 0.28426479,
        0.28407438, 0.28388212, 0.28368811, 0.28349248, 0.28329534,
        0.28309682, 0.28289702, 0.28269608, 0.28249411, 0.28229123,
        0.28208757, 0.28188325, 0.28167839, 0.28147312, 0.28126755,
        0.28106182, 0.28085605, 0.28065036, 0.28044487, 0.28023972,
        0.28003501, 0.27983088, 0.27962745, 0.27942484, 0.27922317,
        0.27902257, 0.27882315, 0.27862504, 0.27842835, 0.27823321,
        0.27803973, 0.27784802, 0.27765821, 0.27747041, 0.27728472,
        0.27710127, 0.27692016, 0.27674151, 0.27656541, 0.27639198,
        0.27622132, 0.27605353, 0.27588871, 0.27572697, 0.2755684 ,
        0.27541309, 0.27526115, 0.27511265, 0.2749677 , 0.27482637,
        0.27468876, 0.27455495, 0.27442501, 0.27429902, 0.27417707,
        0.27405922, 0.27394554, 0.27383611, 0.27373099, 0.27363024,
        0.27353392, 0.2734421 , 0.27335482, 0.27327214, 0.27319411,
        0.27312078, 0.27305219, 0.27298838, 0.27292939, 0.27287526,
        0.27282601, 0.27278168, 0.2727423 , 0.27270789, 0.27267846,
        0.27265404, 0.27263464, 0.27262028, 0.27261096, 0.27260668,
        0.27260746, 0.27261328, 0.27262415, 0.27264006, 0.27266101,
        0.27268696, 0.27271792, 0.27275387, 0.27279477, 0.27284061,
        0.27289136, 0.27294699, 0.27300746, 0.27307275, 0.2731428 ,
        0.27321758]),
 array([0.27437639, 0.27445824, 0.2745447 , 0.2746357 , 0.27473119,
        0.27483112, 0.27493542, 0.27504403, 0.27515689, 0.27527392,
        0.27539506, 0.27552024, 0.27564937, 0.27578239, 0.27591921,
        0.27605974, 0.27620391, 0.27635162, 0.2765028 , 0.27665734,
        0.27681515, 0.27697615, 0.27714022, 0.27730728, 0.27747722,
        0.27764994, 0.27782533, 0.2780033 , 0.27818373, 0.27836651,
        0.27855154, 0.2787387 , 0.27892788, 0.27911897, 0.27931185,
        0.2795064 , 0.27970251, 0.27990006, 0.28009893, 0.280299  ,
        0.28050016, 0.28070227, 0.28090522, 0.28110888, 0.28131313,
        0.28151786, 0.28172293, 0.28192822, 0.28213362, 0.28233898,
        0.2825442 , 0.28274915, 0.2829537 , 0.28315773, 0.28336113,
        0.28356375, 0.28376549, 0.28396622, 0.28416583, 0.28436418,
        0.28456117, 0.28475667, 0.28495056, 0.28514273, 0.28533307,
        0.28552145, 0.28570777, 0.28589191, 0.28607376, 0.28625321,
        0.28643016, 0.28660449, 0.28677611, 0.2869449 , 0.28711077,
        0.28727362, 0.28743335, 0.28758986, 0.28774306, 0.28789285,
        0.28803915, 0.28818187, 0.28832092, 0.28845622, 0.28858768,
        0.28871523, 0.2888388 , 0.28895829, 0.28907366, 0.28918482,
        0.2892917 , 0.28939425, 0.2894924 , 0.28958609, 0.28967527,
        0.28975988, 0.28983987, 0.2899152 , 0.28998581, 0.29005167,
        0.29011273, 0.29016897, 0.29022033, 0.29026681, 0.29030835,
        0.29034495, 0.29037658, 0.29040321, 0.29042484, 0.29044146,
        0.29045304, 0.29045958, 0.29046109, 0.29045755, 0.29044898,
        0.29043537, 0.29041674, 0.29039309, 0.29036444, 0.29033081,
        0.29029221, 0.29024868, 0.29020023, 0.2901469 , 0.29008872,
        0.29002572, 0.28995794, 0.28988543, 0.28980822, 0.28972636,
        0.2896399 , 0.2895489 , 0.2894534 , 0.28935347, 0.28924917,
        0.28914056, 0.2890277 , 0.28891066, 0.28878952, 0.28866434,
        0.2885352 , 0.28840218, 0.28826536, 0.28812483, 0.28798066,
        0.28783294, 0.28768176, 0.28752722, 0.2873694 , 0.28720841,
        0.28704433, 0.28687727, 0.28670733, 0.2865346 , 0.2863592 ,
        0.28618124, 0.2860008 , 0.28581802, 0.28563299, 0.28544583,
        0.28525664, 0.28506555, 0.28487267, 0.28467812, 0.28448201,
        0.28428446, 0.28408558, 0.28388551, 0.28368436, 0.28348225,
        0.2832793 , 0.28307563, 0.28287138, 0.28266665, 0.28246158,
        0.28225629, 0.28205089, 0.28184553, 0.28164031, 0.28143536,
        0.28123081, 0.28102678, 0.28082339, 0.28062076, 0.28041902,
        0.28021829, 0.28001869, 0.27982034, 0.27962335, 0.27942785,
        0.27923396, 0.27904179, 0.27885146, 0.27866308, 0.27847676,
        0.27829262, 0.27811077, 0.27793132, 0.27775438, 0.27758005,
        0.27740844, 0.27723965, 0.27707378, 0.27691093, 0.27675121,
        0.2765947 , 0.2764415 , 0.27629171, 0.27614542, 0.2760027 ,
        0.27586365, 0.27572836, 0.2755969 , 0.27546935, 0.27534579,
        0.27522629, 0.27511093, 0.27499977, 0.27489289, 0.27479034,
        0.2746922 , 0.27459851, 0.27450933, 0.27442472, 0.27434473,
        0.27426941, 0.2741988 , 0.27413294, 0.27407188, 0.27401565,
        0.27396428, 0.27391781, 0.27387626, 0.27383967, 0.27380804,
        0.2737814 , 0.27375977, 0.27374316, 0.27373158, 0.27372504,
        0.27372353, 0.27372707, 0.27373564, 0.27374925, 0.27376788,
        0.27379153, 0.27382018, 0.27385381, 0.2738924 , 0.27393593,
        0.27398438, 0.27403771, 0.27409589, 0.27415889, 0.27422666,
        0.27429918])]
In [135]:
a = 1.0
b = 1.0
c = 0.0
d = 0.0
print("sdeds")
testAB_PS(a, b, c, d, 63, 1000)
testAB_PS(a, b, c, d, 127, 1000)
testAB_PS(a, b, c, d, 255, 1000)
sdeds
(64,)
(64,)
(128,)
(128,)
(256,)
(256,)
Out[135]:
[array([0.28631191, 0.28616154, 0.28600805, 0.28585155, 0.28569211,
        0.28552985, 0.28536484, 0.28519721, 0.28502704, 0.28485444,
        0.28467951, 0.28450236, 0.28432309, 0.28414182, 0.28395865,
        0.28377369, 0.28358706, 0.28339886, 0.28320921, 0.28301822,
        0.28282601, 0.2826327 , 0.2824384 , 0.28224322, 0.28204729,
        0.28185073, 0.28165364, 0.28145616, 0.2812584 , 0.28106047,
        0.2808625 , 0.28066461, 0.28046692, 0.28026954, 0.2800726 ,
        0.27987621, 0.27968049, 0.27948556, 0.27929154, 0.27909854,
        0.27890668, 0.27871607, 0.27852684, 0.27833909, 0.27815293,
        0.27796849, 0.27778587, 0.27760517, 0.27742652, 0.27725001,
        0.27707576, 0.27690386, 0.27673443, 0.27656756, 0.27640335,
        0.27624191, 0.27608333, 0.2759277 , 0.27577512, 0.27562569,
        0.27547948, 0.2753366 , 0.27519712, 0.27506113, 0.27492872,
        0.27479995, 0.27467491, 0.27455368, 0.27443633, 0.27432292,
        0.27421353, 0.27410822, 0.27400706, 0.2739101 , 0.27381741,
        0.27372904, 0.27364505, 0.27356547, 0.27349038, 0.27341979,
        0.27335377, 0.27329235, 0.27323557, 0.27318346, 0.27313605,
        0.27309337, 0.27305545, 0.27302231, 0.27299396, 0.27297044,
        0.27295174, 0.27293789, 0.27292889, 0.27292475, 0.27292546,
        0.27293103, 0.27294146, 0.27295673, 0.27297684, 0.27300178,
        0.27303154, 0.27306608, 0.2731054 , 0.27314947, 0.27319827,
        0.27325175, 0.2733099 , 0.27337268, 0.27344004, 0.27351195,
        0.27358837, 0.27366924, 0.27375453, 0.27384417, 0.27393812,
        0.27403632, 0.2741387 , 0.27424522, 0.27435579, 0.27447037,
        0.27458887, 0.27471122, 0.27483736, 0.27496721, 0.27510068,
        0.2752377 , 0.27537819, 0.27552206, 0.27566922, 0.27581958,
        0.27597306, 0.27612957, 0.276289  , 0.27645126, 0.27661626,
        0.27678389, 0.27695406, 0.27712666, 0.27730159, 0.27747873,
        0.277658  , 0.27783927, 0.27802244, 0.2782074 , 0.27839403,
        0.27858223, 0.27877188, 0.27896286, 0.27915507, 0.27934838,
        0.27954268, 0.27973785, 0.27993378, 0.28013035, 0.28032743,
        0.28052491, 0.28072268, 0.2809206 , 0.28111857, 0.28131646,
        0.28151415, 0.28171153, 0.28190847, 0.28210486, 0.28230058,
        0.28249551, 0.28268954, 0.28288254, 0.2830744 , 0.283265  ,
        0.28345424, 0.28364199, 0.28382815, 0.28401259, 0.28419522,
        0.28437591, 0.28455457, 0.28473108, 0.28490533, 0.28507723,
        0.28524667, 0.28541354, 0.28557775, 0.2857392 , 0.28589778,
        0.28605341, 0.28620599, 0.28635543, 0.28650163, 0.28664452,
        0.286784  , 0.28691999, 0.28705241, 0.28718118, 0.28730622,
        0.28742745, 0.28754481, 0.28765822, 0.28776761, 0.28787292,
        0.28797409, 0.28807104, 0.28816374, 0.28825211, 0.28833611,
        0.28841568, 0.28849078, 0.28856137, 0.28862739, 0.28868881,
        0.2887456 , 0.28879771, 0.28884512, 0.2888878 , 0.28892572,
        0.28895887, 0.28898721, 0.28901073, 0.28902943, 0.28904328,
        0.28905228, 0.28905643, 0.28905572, 0.28905015, 0.28903972,
        0.28902445, 0.28900433, 0.28897939, 0.28894964, 0.28891509,
        0.28887577, 0.2888317 , 0.2887829 , 0.28872942, 0.28867127,
        0.28860849, 0.28854112, 0.28846921, 0.28839279, 0.28831192,
        0.28822663, 0.28813698, 0.28804303, 0.28794483, 0.28784245,
        0.28773593, 0.28762535, 0.28751078, 0.28739228, 0.28726992,
        0.28714378, 0.28701393, 0.28688045, 0.28674343, 0.28660294,
        0.28645907]),
 array([0.28657219, 0.28641911, 0.28626277, 0.28610325, 0.28594065,
        0.28577507, 0.2856066 , 0.28543536, 0.28526144, 0.28508494,
        0.28490598, 0.28472466, 0.28454109, 0.28435539, 0.28416765,
        0.28397801, 0.28378656, 0.28359343, 0.28339873, 0.28320258,
        0.2830051 , 0.2828064 , 0.28260662, 0.28240585, 0.28220424,
        0.28200189, 0.28179894, 0.2815955 , 0.28139169, 0.28118764,
        0.28098348, 0.28077931, 0.28057528, 0.28037149, 0.28016808,
        0.27996517, 0.27976287, 0.27956131, 0.27936061, 0.2791609 ,
        0.27896228, 0.27876489, 0.27856884, 0.27837425, 0.27818123,
        0.27798991, 0.27780039, 0.27761279, 0.27742723, 0.27724382,
        0.27706266, 0.27688387, 0.27670756, 0.27653382, 0.27636277,
        0.2761945 , 0.27602912, 0.27586674, 0.27570744, 0.27555132,
        0.27539848, 0.275249  , 0.27510299, 0.27496052, 0.27482168,
        0.27468656, 0.27455524, 0.27442779, 0.2743043 , 0.27418483,
        0.27406946, 0.27395827, 0.27385131, 0.27374864, 0.27365035,
        0.27355647, 0.27346707, 0.2733822 , 0.27330192, 0.27322626,
        0.27315529, 0.27308903, 0.27302753, 0.27297083, 0.27291896,
        0.27287195, 0.27282983, 0.27279263, 0.27276037, 0.27273306,
        0.27271073, 0.27269338, 0.27268104, 0.2726737 , 0.27267136,
        0.27267404, 0.27268173, 0.27269442, 0.27271211, 0.27273479,
        0.27276244, 0.27279505, 0.27283259, 0.27287505, 0.27292239,
        0.2729746 , 0.27303163, 0.27309346, 0.27316004, 0.27323134,
        0.27330732, 0.27338792, 0.27347311, 0.27356282, 0.273657  ,
        0.2737556 , 0.27385857, 0.27396582, 0.27407731, 0.27419296,
        0.27431271, 0.27443648, 0.2745642 , 0.27469578, 0.27483116,
        0.27497025, 0.27511297, 0.27525923, 0.27540894, 0.27556201,
        0.27571835, 0.27587787, 0.27604046, 0.27620604, 0.2763745 ,
        0.27654575, 0.27671966, 0.27689616, 0.27707511, 0.27725643,
        0.27744   , 0.2776257 , 0.27781343, 0.27800308, 0.27819452,
        0.27838765, 0.27858235, 0.27877849, 0.27897597, 0.27917467,
        0.27937445, 0.27957522, 0.27977683, 0.27997917, 0.28018213,
        0.28038557, 0.28058937, 0.28079342, 0.28099759, 0.28120175,
        0.28140578, 0.28160957, 0.28181298, 0.2820159 , 0.2822182 ,
        0.28241976, 0.28262045, 0.28282017, 0.28301879, 0.28321618,
        0.28341223, 0.28360683, 0.28379984, 0.28399117, 0.28418069,
        0.28436829, 0.28455385, 0.28473726, 0.28491842, 0.28509722,
        0.28527354, 0.28544728, 0.28561833, 0.2857866 , 0.28595198,
        0.28611437, 0.28627367, 0.28642979, 0.28658264, 0.28673211,
        0.28687813, 0.2870206 , 0.28715944, 0.28729457, 0.28742589,
        0.28755334, 0.28767684, 0.28779631, 0.28791168, 0.28802288,
        0.28812984, 0.2882325 , 0.28833081, 0.28842469, 0.28851409,
        0.28859896, 0.28867925, 0.2887549 , 0.28882588, 0.28889214,
        0.28895364, 0.28901034, 0.28906221, 0.28910922, 0.28915134,
        0.28918854, 0.28922081, 0.28924812, 0.28927045, 0.2892878 ,
        0.28930014, 0.28930749, 0.28930982, 0.28930714, 0.28929945,
        0.28928676, 0.28926907, 0.28924639, 0.28921874, 0.28918613,
        0.28914859, 0.28910613, 0.28905879, 0.28900658, 0.28894954,
        0.28888772, 0.28882113, 0.28874983, 0.28867385, 0.28859324,
        0.28850806, 0.28841835, 0.28832416, 0.28822555, 0.28812259,
        0.28801533, 0.28790384, 0.28778819, 0.28766844, 0.28754466,
        0.28741694, 0.28728535, 0.28714997, 0.28701088, 0.28686816,
        0.2867219 ]),
 array([0.28771039, 0.28755647, 0.28739927, 0.28723886, 0.28707535,
        0.28690884, 0.28673943, 0.28656723, 0.28639232, 0.28621483,
        0.28603485, 0.2858525 , 0.28566788, 0.28548111, 0.2852923 ,
        0.28510156, 0.284909  , 0.28471476, 0.28451893, 0.28432164,
        0.284123  , 0.28392315, 0.28372219, 0.28352025, 0.28331744,
        0.28311391, 0.28290975, 0.2827051 , 0.28250009, 0.28229483,
        0.28208944, 0.28188406, 0.28167881, 0.2814738 , 0.28126916,
        0.28106502, 0.28086151, 0.28065873, 0.28045681, 0.28025588,
        0.28005606, 0.27985746, 0.27966021, 0.27946443, 0.27927022,
        0.27907772, 0.27888703, 0.27869828, 0.27851157, 0.27832701,
        0.27814472, 0.27796482, 0.27778739, 0.27761256, 0.27744043,
        0.2772711 , 0.27710468, 0.27694126, 0.27678094, 0.27662382,
        0.27647   , 0.27631956, 0.2761726 , 0.27602921, 0.27588947,
        0.27575346, 0.27562127, 0.27549298, 0.27536867, 0.2752484 ,
        0.27513226, 0.27502031, 0.27491262, 0.27480926, 0.27471028,
        0.27461575, 0.27452572, 0.27444025, 0.27435939, 0.27428319,
        0.27421169, 0.27414493, 0.27408297, 0.27402583, 0.27397355,
        0.27392616, 0.27388368, 0.27384616, 0.2738136 , 0.27378602,
        0.27376345, 0.2737459 , 0.27373337, 0.27372588, 0.27372343,
        0.27372602, 0.27373365, 0.27374631, 0.273764  , 0.27378671,
        0.27381442, 0.27384712, 0.27388478, 0.27392739, 0.27397491,
        0.27402733, 0.2740846 , 0.2741467 , 0.27421358, 0.27428521,
        0.27436154, 0.27444253, 0.27452812, 0.27461828, 0.27471293,
        0.27481203, 0.27491551, 0.27502332, 0.27513539, 0.27525164,
        0.27537202, 0.27549644, 0.27562484, 0.27575713, 0.27589324,
        0.27603309, 0.27617658, 0.27632364, 0.27647417, 0.27662808,
        0.27678529, 0.27694569, 0.2771092 , 0.2772757 , 0.27744511,
        0.27761732, 0.27779222, 0.27796971, 0.27814968, 0.27833203,
        0.27851665, 0.27870342, 0.27889223, 0.27908297, 0.27927552,
        0.27946977, 0.27966559, 0.27986288, 0.28006151, 0.28026137,
        0.28046233, 0.28066427, 0.28086707, 0.28107061, 0.28127476,
        0.28147941, 0.28168442, 0.28188968, 0.28209507, 0.28230045,
        0.28250571, 0.28271071, 0.28291535, 0.28311949, 0.28332301,
        0.28352579, 0.2837277 , 0.28392863, 0.28412846, 0.28432706,
        0.28452431, 0.2847201 , 0.2849143 , 0.2851068 , 0.28529749,
        0.28548625, 0.28567296, 0.28585752, 0.28603981, 0.28621972,
        0.28639715, 0.28657198, 0.28674411, 0.28691344, 0.28707987,
        0.28724329, 0.28740361, 0.28756073, 0.28771456, 0.287865  ,
        0.28801196, 0.28815536, 0.28829511, 0.28843111, 0.2885633 ,
        0.2886916 , 0.28881591, 0.28893618, 0.28905233, 0.28916428,
        0.28927197, 0.28937534, 0.28947432, 0.28956885, 0.28965888,
        0.28974435, 0.28982522, 0.28990142, 0.28997292, 0.29003967,
        0.29010164, 0.29015878, 0.29021107, 0.29025846, 0.29030093,
        0.29033846, 0.29037102, 0.2903986 , 0.29042117, 0.29043872,
        0.29045125, 0.29045874, 0.29046119, 0.2904586 , 0.29045097,
        0.29043831, 0.29042062, 0.29039791, 0.2903702 , 0.2903375 ,
        0.29029983, 0.29025723, 0.2902097 , 0.29015728, 0.29010001,
        0.29003791, 0.28997103, 0.2898994 , 0.28982306, 0.28974207,
        0.28965648, 0.28956632, 0.28947167, 0.28937256, 0.28926908,
        0.28916127, 0.2890492 , 0.28893294, 0.28881256, 0.28868814,
        0.28855974, 0.28842744, 0.28829133, 0.28815148, 0.28800799,
        0.28786093])]
In [136]:
a = 1.0
b = 0.01
c = 0.0; 
d = 0.0
print("ssfw")
testAB_PS(a, b, c, d, 127, 1000)
testAB_PS(a, b, c, d, 127, 10000)
ssfw
(128,)
(128,)
(128,)
(128,)
Out[136]:
[array([4.73207819e-01, 4.37737562e-01, 4.03181822e-01, 3.69751986e-01,
        3.37628998e-01, 3.06962852e-01, 2.77872745e-01, 2.50447854e-01,
        2.24748650e-01, 2.00808665e-01, 1.78636642e-01, 1.58218964e-01,
        1.39522280e-01, 1.22496241e-01, 1.07076276e-01, 9.31863210e-02,
        8.07414596e-02, 6.96504053e-02, 5.98177967e-02, 5.11462708e-02,
        4.35382942e-02, 3.68977425e-02, 3.11312252e-02, 2.61491609e-02,
        2.18666156e-02, 1.82039183e-02, 1.50870761e-02, 1.24480099e-02,
        1.02246347e-02, 8.36080878e-03, 6.80617670e-03, 5.51592684e-03,
        4.45048654e-03, 3.57517306e-03, 2.85981769e-03, 2.27837773e-03,
        1.80854860e-03, 1.43138628e-03, 1.13094819e-03, 8.93958455e-04,
        7.09502197e-04, 5.68751831e-04, 4.64727408e-04, 3.92092153e-04,
        3.46983810e-04, 3.26882102e-04, 3.30512466e-04, 3.57786210e-04,
        4.09777300e-04, 4.88735962e-04, 5.98139245e-04, 7.42778492e-04,
        9.28883270e-04, 1.16428089e-03, 1.45858998e-03, 1.82344599e-03,
        2.27275577e-03, 2.82297780e-03, 3.49342381e-03, 4.30657718e-03,
        5.28842215e-03, 6.46877746e-03, 7.88162606e-03, 9.56543144e-03,
        1.15634289e-02, 1.39238783e-02, 1.67002617e-02, 1.99514098e-02,
        2.37415359e-02, 2.81401583e-02, 3.32218902e-02, 3.90660760e-02,
        4.57562539e-02, 5.33794287e-02, 6.20251383e-02, 7.17843024e-02,
        8.27478486e-02, 9.50051125e-02, 1.08642019e-01, 1.23739061e-01,
        1.40369086e-01, 1.58594942e-01, 1.78467002e-01, 2.00020627e-01,
        2.23273623e-01, 2.48223757e-01, 2.74846411e-01, 3.03092437e-01,
        3.32886306e-01, 3.64124629e-01, 3.96675119e-01, 4.30376084e-01,
        4.65036495e-01, 5.00436702e-01, 5.36329820e-01, 5.72443817e-01,
        6.08484305e-01, 6.44138011e-01, 6.79076891e-01, 7.12962831e-01,
        7.45452835e-01, 7.76204614e-01, 8.04882445e-01, 8.31163168e-01,
        8.54742166e-01, 8.75339190e-01, 8.92703866e-01, 9.06620729e-01,
        9.16913656e-01, 9.23449562e-01, 9.26141251e-01, 9.24949347e-01,
        9.19883225e-01, 9.11000939e-01, 8.98408119e-01, 8.82255888e-01,
        8.62737835e-01, 8.40086153e-01, 8.14567028e-01, 7.86475417e-01,
        7.56129364e-01, 7.23863994e-01, 6.90025357e-01, 6.54964273e-01,
        6.19030328e-01, 5.82566171e-01, 5.45902236e-01, 5.09351993e-01]),
 array([4.92427555e-01, 4.56803327e-01, 4.21976181e-01, 3.88166776e-01,
        3.55566225e-01, 3.24335387e-01, 2.94604813e-01, 2.66475315e-01,
        2.40019070e-01, 2.15281200e-01, 1.92281735e-01, 1.71017874e-01,
        1.51466466e-01, 1.33586624e-01, 1.17322397e-01, 1.02605426e-01,
        8.93575285e-02, 7.74931588e-02, 6.69216981e-02, 5.75495442e-02,
        4.92819799e-02, 4.20248049e-02, 3.56857263e-02, 3.01755108e-02,
        2.54089056e-02, 2.13053425e-02, 1.77894422e-02, 1.47913380e-02,
        1.22468416e-02, 1.00974729e-02, 8.29037606e-03, 6.77814400e-03,
        5.51857197e-03, 4.47435837e-03, 3.61277044e-03, 2.90528944e-03,
        2.32724808e-03, 1.85747057e-03, 1.47792389e-03, 1.17338676e-03,
        9.31141292e-04, 7.40691032e-04, 5.93508025e-04, 4.82811078e-04,
        4.03376842e-04, 3.51385238e-04, 3.24300538e-04, 3.20789286e-04,
        3.40675938e-04, 3.84936759e-04, 4.55732063e-04, 5.56476398e-04,
        6.91945922e-04, 8.68421848e-04, 1.09386867e-03, 1.37814567e-03,
        1.73325015e-03, 2.17359031e-03, 2.71628535e-03, 3.38148950e-03,
        4.19273535e-03, 5.17729055e-03, 6.36652005e-03, 7.79624401e-03,
        9.50707956e-03, 1.15447520e-02, 1.39603596e-02, 1.68105729e-02,
        2.01577494e-02, 2.40699423e-02, 2.86207791e-02, 3.38891912e-02,
        3.99589682e-02, 4.69181207e-02, 5.48580291e-02, 6.38723676e-02,
        7.40557908e-02, 8.55023790e-02, 9.83038468e-02, 1.12547522e-01,
        1.28314119e-01, 1.45675331e-01, 1.64691281e-01, 1.85407887e-01,
        2.07854186e-01, 2.32039699e-01, 2.57951900e-01, 2.85553877e-01,
        3.14782266e-01, 3.45545541e-01, 3.77722747e-01, 4.11162757e-01,
        4.45684118e-01, 4.81075542e-01, 5.17097098e-01, 5.53482127e-01,
        5.89939875e-01, 6.26158855e-01, 6.61810873e-01, 6.96555683e-01,
        7.30046167e-01, 7.61933956e-01, 7.91875347e-01, 8.19537406e-01,
        8.44604084e-01, 8.66782202e-01, 8.85807145e-01, 9.01448109e-01,
        9.13512760e-01, 9.21851169e-01, 9.26358913e-01, 9.26979250e-01,
        9.23704310e-01, 9.16575258e-01, 9.05681437e-01, 8.91158501e-01,
        8.73185604e-01, 8.51981725e-01, 8.27801226e-01, 8.00928779e-01,
        7.71673800e-01, 7.40364537e-01, 7.07341977e-01, 6.72953726e-01,
        6.37548006e-01, 6.01467931e-01, 5.65046164e-01, 5.28600086e-01]),
 array([5.15256003e-01, 4.79379122e-01, 4.44147266e-01, 4.09795035e-01,
        3.76528812e-01, 3.44525646e-01, 3.13932819e-01, 2.84868026e-01,
        2.57420135e-01, 2.31650447e-01, 2.07594381e-01, 1.85263501e-01,
        1.64647811e-01, 1.45718226e-01, 1.28429152e-01, 1.12721088e-01,
        9.85232103e-02, 8.57558503e-02, 7.43328495e-02, 6.41637337e-02,
        5.51556866e-02, 4.72153035e-02, 4.02501124e-02, 3.41698628e-02,
        2.88875825e-02, 2.43204162e-02, 2.03902554e-02, 1.70241818e-02,
        1.41547394e-02, 1.17200602e-02, 9.66386347e-03, 7.93534993e-03,
        6.48901281e-03, 5.28438438e-03, 4.28573612e-03, 3.46174847e-03,
        2.78516406e-03, 2.23243605e-03, 1.78338140e-03, 1.42084644e-03,
        1.13039103e-03, 8.99995328e-04, 7.19792903e-04, 5.81832486e-04,
        4.79870752e-04, 4.09197984e-04, 3.66498526e-04, 3.49747708e-04,
        3.58146721e-04, 3.92096419e-04, 4.53210502e-04, 5.44367847e-04,
        6.69803164e-04, 8.35234658e-04, 1.04802704e-03, 1.31738812e-03,
        1.65459692e-03, 2.07326147e-03, 2.58960353e-03, 3.22276748e-03,
        3.99514911e-03, 4.93273905e-03, 6.06547377e-03, 7.42758516e-03,
        9.05793768e-03, 1.10003399e-02, 1.33038153e-02, 1.60228145e-02,
        1.92173510e-02, 2.29530386e-02, 2.73010109e-02, 3.23376989e-02,
        3.81444475e-02, 4.48069485e-02, 5.24144739e-02, 6.10588916e-02,
        7.08334546e-02, 8.18313568e-02, 9.41440555e-02, 1.07859369e-01,
        1.23059365e-01, 1.39818064e-01, 1.58198995e-01, 1.78252644e-01,
        2.00013846e-01, 2.23499195e-01, 2.48704525e-01, 2.75602550e-01,
        3.04140741e-01, 3.34239512e-01, 3.65790817e-01, 3.98657217e-01,
        4.32671495e-01, 4.67636886e-01, 5.03327961e-01, 5.39492208e-01,
        5.75852315e-01, 6.12109161e-01, 6.47945478e-01, 6.83030141e-01,
        7.17023015e-01, 7.49580267e-01, 7.80360028e-01, 8.09028290e-01,
        8.35264882e-01, 8.58769393e-01, 8.79266881e-01, 8.96513212e-01,
        9.10299899e-01, 9.20458298e-01, 9.26863033e-01, 9.29434581e-01,
        9.28140907e-01, 9.22998135e-01, 9.14070207e-01, 9.01467567e-01,
        8.85344881e-01, 8.65897880e-01, 8.43359407e-01, 8.17994781e-01,
        7.90096601e-01, 7.59979147e-01, 7.27972511e-01, 6.94416619e-01,
        6.59655296e-01, 6.24030509e-01, 5.87876924e-01, 5.51516901e-01])]

With First Order Term

In [137]:
def testABC_PS(a, b, c, d, nx, nt):

    L = np.pi
    T = 4.0


    x = np.linspace(-L, L, nx+1)
    t = np.linspace(0, T, nt+1)
    
    u = np.exp(-x**2)#np.sin(np.pi*xin)
    ui = u;
    U_CDI = np.zeros((nx+1, nt+1))
    error_2norm_CDI_PSTD = np.zeros(nt+1)

    # Discretization of space and time
    dx = 2*L/(nx+1)
    # dt = 0.1*dx**2/D#T/(nt-1)
    # dx = L/(nx+1)
    dt = T/(nt+1)


    #plt.figure()
    #plt.plot(u, label='Initial')
    def numerical_scheme_CDI(u_CDI, a, b, c, d, dx, dt):
        alpha = a*dt/dx
        beta = b*dt/dx**2
        nx = len(u_CDI)-1;
        A = np.zeros((nx+1, nx+1))
        bi = d*dt+u_CDI;
        for j in range(1, len(u_CDI)-1):

            A[j, j-1] = -(beta-alpha/2)
            A[j, j] = 1+2*beta-c*dt
            A[j, j+1] = -(beta+alpha/2)

        # Applying first order Periodic Boundary Condition
        A[0, 0]=(1+2*beta);
        A[0, 1]= (-alpha/2-beta)
        A[0, nx]= (alpha/2-beta)    


        A[nx, 0]= (-alpha/2-beta) ;
        A[nx, nx-1]= (alpha/2 - beta)
        A[nx, nx]= 1+2*beta    


        # Applying second order periodic boundary conditions
        # bt = d*dt+u_CDI;
        # kk = np.zeros(2)
        # bi = [*kk, *bt]

        # A[0, 0]=0;
        # A[0, 1]= 1
        # A[0, N]= -1
        # A[0, N+1]=0

        # A[1, 0] = -1
        # A[1,1] = 0
        # A[1, 2] = 1
        # A[1, N-1] = 1
        # A[1, N] = 0 
        # A[1, N+1] = -1


        # Boundary condition
        # A[0, 0] = -3
        # A[0, 1] = 4
        # A[0, 2] = -1
        # bi[0] = 0

        # A[N, N] = 3
        # A[N, N-1] = -4
        # A[N, N-2] = 1
        # bi[N] = 0


        # A[0,0] = 1.0
        # # A[N,N] = 1.0

        # bi[0] = l
        # bi[N] = r


        u_CDI = solve(A, bi)

        return u_CDI



    # Apply the numerical scheme
    # u = np.exp(-((x+0.5)**2)/(0.00125))
    U_CDI = np.zeros((nx+1, nt+1))
    error_2norm_QI_PSTD = np.zeros(nt)

    u_CDI=u;
    plt.plot(x,u)
    for n in range(nt):
        u_CDI = numerical_scheme_CDI(u_CDI, a, b, c, d, dx, dt)
        U_CDI[:, n] = u_CDI

    #plt.plot(u_CDI, label='Final')
    #plt.show()
    #plt.close()

    plt.rcParams['font.size'] = 18
    fig = plt.figure() 
    plt.imshow(U_CDI, cmap='viridis', extent=[0, T, 0, 2*L], aspect='auto')
    plt.colorbar()
    title = "CDI - nx = "+str(nx+1)
    plt.title(title)
    plt.xlabel('Time')
    plt.ylabel('Position')
    plt.show()


    k = 2*np.pi*np.fft.rfftfreq(nx+1, 2*L/(nx))
    k2=k**2;

    # Defining variables
    U_ps = np.zeros((nx+1, nt+1))
    U_psk = np.zeros((nx+1, nt+1))
    #plt.figure()
    #plt.plot(xin, u, color = 'b', label = "initial")
    # left = (0.025/(np.sqrt(0.000625 + 0.02*t)))*(np.exp((-(0.5 - t)**2)/(0.00125 + 0.04*t)))
    # right = (0.025/(np.sqrt(0.000625 + 0.02*t)))*(np.exp((-(1.5 - t)**2)/(0.00125 + 0.04*t)))

    # Solving over time
    for i in range(nt+1): # 
        
        uk = np.fft.rfft(u)

        uk[:] = (uk[:])/(1-(-b*k**2 + a*1j*k+c)*dt) 
        # print(np.shape(uk))    
        u = np.fft.irfft(uk) #+ d*dt
        U_ps[:, i] = u
        error =  U_CDI[:, i] - u
        error_2norm_CDI_PSTD[i] =np.linalg.norm(error)/np.sqrt(nx)


    fig = plt.figure()  
    plt.imshow((U_ps), cmap='viridis',  extent=[0, T, 0, 2*L], aspect='auto')
    plt.colorbar()
    plt.title("CD Scheme with PSTD")
    plt.xlabel('Time')
    plt.ylabel('Position')
    title = "PSTDI - nx = "+str(nx+1)
    plt.title(title)
    plt.grid(True)
    plt.show()
    plt.close()
    
    U_psk = np.zeros((nx+1, nt+1))
    error_2norm_PSTDI_PSTDK = np.zeros(nt+1)
    
    error_2norm_CDI_PSTDK = np.zeros(nt+1)
    # Spatial grid
    m=nx+1                            # Number of grid points in space
    L = np.pi                   # Width of spatial domain
    x = np.arange(-m/2,m/2)*(2*L/m)   # Grid points
    dx = x[1]-x[0]                  # Grid spacing

    # Temporal grid
    tmax=T     # Final time
    N = nt+1       # number grid points in time
    dt = tmax/N   # interval between output times

    xi = np.fft.fftfreq(m)*m*np.pi/L  # Wavenumber "grid"
    xi = 2*np.pi*np.fft.fftfreq(nx+1, 2*L/(nx))
    
    # k = 2*np.pi*np.fft.rfftfreq(nx+1, 2*L/(nx))
    
    
    # (this is the order in which numpy's FFT gives the frequencies)

    # Initial data
    u = np.exp(-x**2)# np.sin(2*x)**2 * (x<-L/4)
    uhat0 = np.fft.fft(u)
    plt.plot(x,u)
    epsilon=b  # Diffusion coefficient
    # a = 1.0       # Advection coefficient

    print(np.shape(xi))
    print(np.shape(u))
    # Now we solve the problem
    for i in range(0,N):
        t = i*dt
        uhat = np.exp((1.j*xi*a - epsilon*xi**2+c)*t) * uhat0
        u = np.real(np.fft.ifft(uhat)) #+ d*t
        
        U_psk[:, i] = u
        error =  U_CDI[:, i] - u
        error_2norm_CDI_PSTDK[i] =np.linalg.norm(error)/np.sqrt(nx)
        error =  U_ps[:, i] - u
        error_2norm_PSTDI_PSTDK[i] =np.linalg.norm(error)/np.sqrt(nx)

    
    fig = plt.figure()  
    plt.imshow((U_psk), cmap='viridis',  extent=[0, T, 0, 2*L], aspect='auto')
    plt.colorbar()
    plt.title("CD Scheme with PSTDK")
    plt.xlabel('Time')
    plt.ylabel('Position')
    title = "PSTDK - nx = "+str(nx+1)
    plt.title(title)
    plt.grid(True)
    plt.show()
    plt.close()
    
    plt.figure();
    # plt.plot(x, ui, '-', color='g', label = 'Initial'); 
    plt.plot(x, U_CDI[:, int(nt/2)], '-', color='b', label = 'CDI Method'); 
    plt.plot(x, U_ps[:, int(nt/2)], '-', color = 'r', label = 'Spectral Method');
    plt.plot(x, U_psk[:, int(nt/2)+1], '-', color = 'g', label = 'Ketchinson Spectral Method');
    
    res =  [U_CDI[:, nt-1], U_ps[:, nt-1], U_psk[:, nt]]
    
    plt.legend(loc="upper right")
    title = "Comparison - nx = "+str(nx+1)
    plt.title(title)
    plt.show()
    plt.close()
    plt.figure()
    plt.plot(error_2norm_CDI_PSTD[0:nt], '.', label = 'PSTDI_CDI')
    
    plt.plot(error_2norm_CDI_PSTDK[0:nt], '.', label = 'PSTDK_CDI')
    
    plt.plot(error_2norm_PSTDI_PSTDK[0:nt], '.', label = 'PSTDI_PSTDK')
    plt.legend()
    plt.show()
    plt.close()
    return res
In [138]:
a = 1.0
b = 0.01
c = 0.5; 
d = 0.0

testABC_PS(a, b, c, d, 63, 1000)
testABC_PS(a, b, c, d, 127, 1000)
res = testABC_PS(a, b, c, d, 255, 1000)

plt.figure();
plt.plot(res[0], '-', color='b', label = 'CDI Method'); 
plt.plot(res[1], '-', color = 'r', label = 'Spectral Method');
plt.plot(res[2], '-', color = 'g', label = 'Ketchinson Spectral Method');
plt.legend(loc="upper right")
title = "Comparison - nx = "+str(255+1)
plt.title(title)
plt.show()
plt.close()
(64,)
(64,)
(128,)
(128,)
(256,)
(256,)
In [139]:
a = 1.0
b = 1.0
c = 0.5; 
d = 0.0

testABC_PS(a, b, c, d, 63, 1000)
testABC_PS(a, b, c, d, 127, 1000)
res = testABC_PS(a, b, c, d, 255, 1000)

plt.figure();
plt.plot(res[0], '-', color='b', label = 'CDI Method'); 
plt.plot(res[1], '-', color = 'r', label = 'Spectral Method');
plt.plot(res[2], '-', color = 'g', label = 'Ketchinson Spectral Method');
plt.legend(loc="upper right")
title = "Comparison - nx = "+str(255+1)
plt.title(title)
plt.show()
plt.close()
(64,)
(64,)
(128,)
(128,)
(256,)
(256,)
In [140]:
a = 0.1
b = 1.0
c = 0.5; 
d = 0.0

testABC_PS(a, b, c, d, 63, 1000)
testABC_PS(a, b, c, d, 127, 1000)
res = testABC_PS(a, b, c, d, 255, 1000)

plt.figure();
plt.plot(res[0], '-', color='b', label = 'CDI Method'); 
plt.plot(res[1], '-', color = 'r', label = 'Spectral Method');
plt.plot(res[2], '-', color = 'g', label = 'Ketchinson Spectral Method');
plt.legend(loc="upper right")
title = "Comparison - nx = "+str(255+1)
plt.title(title)
plt.show()
plt.close()
(64,)
(64,)
(128,)
(128,)
(256,)
(256,)

With Zeroth and First Order Terms

In [141]:
def test_PS(a, b, c, d, nx, nt):

    L = np.pi
    T = 4.0


    x = np.linspace(-L, L, nx+1)
    t = np.linspace(0, T, nt+1)
    
    u = np.exp(-x**2)#np.sin(np.pi*xin)
    ui = u;
    U_CDI = np.zeros((nx+1, nt+1))
    error_2norm_CDI_PSTD = np.zeros(nt+1)

    # Discretization of space and time
    dx = 2*L/(nx+1)
    # dt = 0.1*dx**2/D#T/(nt-1)
    # dx = L/(nx+1)
    dt = T/(nt+1)


    #plt.figure()
    #plt.plot(u, label='Initial')
    def numerical_scheme_CDI(u_CDI, a, b, c, d, dx, dt):
        alpha = a*dt/dx
        beta = b*dt/dx**2
        nx = len(u_CDI)-1;
        A = np.zeros((nx+1, nx+1))
        bi = d*dt+u_CDI;
        for j in range(1, len(u_CDI)-1):

            A[j, j-1] = -(beta-alpha/2)
            A[j, j] = 1+2*beta-c*dt
            A[j, j+1] = -(beta+alpha/2)

        # Applying first order Periodic Boundary Condition
        A[0, 0]=(1+2*beta);
        A[0, 1]= (-alpha/2-beta)
        A[0, nx]= (alpha/2-beta)    


        A[nx, 0]= (-alpha/2-beta) ;
        A[nx, nx-1]= (alpha/2 - beta)
        A[nx, nx]= 1+2*beta    


        # Applying second order periodic boundary conditions
        # bt = d*dt+u_CDI;
        # kk = np.zeros(2)
        # bi = [*kk, *bt]

        # A[0, 0]=0;
        # A[0, 1]= 1
        # A[0, N]= -1
        # A[0, N+1]=0

        # A[1, 0] = -1
        # A[1,1] = 0
        # A[1, 2] = 1
        # A[1, N-1] = 1
        # A[1, N] = 0 
        # A[1, N+1] = -1


        # Boundary condition
        # A[0, 0] = -3
        # A[0, 1] = 4
        # A[0, 2] = -1
        # bi[0] = 0

        # A[N, N] = 3
        # A[N, N-1] = -4
        # A[N, N-2] = 1
        # bi[N] = 0


        # A[0,0] = 1.0
        # # A[N,N] = 1.0

        # bi[0] = l
        # bi[N] = r


        u_CDI = solve(A, bi)

        return u_CDI



    # Apply the numerical scheme
    # u = np.exp(-((x+0.5)**2)/(0.00125))
    U_CDI = np.zeros((nx+1, nt+1))
    error_2norm_QI_PSTD = np.zeros(nt)

    u_CDI=u;
    plt.plot(x,u)
    for n in range(nt):
        u_CDI = numerical_scheme_CDI(u_CDI, a, b, c, d, dx, dt)
        U_CDI[:, n] = u_CDI

    #plt.plot(u_CDI, label='Final')
    #plt.show()
    #plt.close()

    plt.rcParams['font.size'] = 18
    fig = plt.figure() 
    plt.imshow(U_CDI, cmap='viridis', extent=[0, T, 0, 2*L], aspect='auto')
    plt.colorbar()
    title = "CDI - nx = "+str(nx+1)
    plt.title(title)
    plt.xlabel('Time')
    plt.ylabel('Position')
    plt.show()


    k = 2*np.pi*np.fft.rfftfreq(nx+1, 2*L/(nx))
    k2=k**2;

    # Defining variables
    U_ps = np.zeros((nx+1, nt+1))
    U_psk = np.zeros((nx+1, nt+1))
    #plt.figure()
    #plt.plot(xin, u, color = 'b', label = "initial")
    # left = (0.025/(np.sqrt(0.000625 + 0.02*t)))*(np.exp((-(0.5 - t)**2)/(0.00125 + 0.04*t)))
    # right = (0.025/(np.sqrt(0.000625 + 0.02*t)))*(np.exp((-(1.5 - t)**2)/(0.00125 + 0.04*t)))

    # Solving over time
    for i in range(nt+1): # 
        
        uk = np.fft.rfft(u)

        uk[:] = (uk[:])/(1-(-b*k**2 + a*1j*k+c)*dt) 
        # print(np.shape(uk))    
        u = np.fft.irfft(uk) + d*dt
        U_ps[:, i] = u
        error =  U_CDI[:, i] - u
        error_2norm_CDI_PSTD[i] =np.linalg.norm(error)/np.sqrt(nx)


    fig = plt.figure()  
    plt.imshow((U_ps), cmap='viridis',  extent=[0, T, 0, 2*L], aspect='auto')
    plt.colorbar()
    plt.title("CD Scheme with PSTD")
    plt.xlabel('Time')
    plt.ylabel('Position')
    title = "PSTDI - nx = "+str(nx+1)
    plt.title(title)
    plt.grid(True)
    plt.show()
    plt.close()
    
    U_psk = np.zeros((nx+1, nt+1))
    error_2norm_PSTDI_PSTDK = np.zeros(nt+1)
    
    error_2norm_CDI_PSTDK = np.zeros(nt+1)
    # Spatial grid
    m=nx+1                            # Number of grid points in space
    L = np.pi                   # Width of spatial domain
    x = np.arange(-m/2,m/2)*(2*L/m)   # Grid points
    dx = x[1]-x[0]                  # Grid spacing

    # Temporal grid
    tmax=T     # Final time
    N = nt+1       # number grid points in time
    dt = tmax/N   # interval between output times

    xi = np.fft.fftfreq(m)*m*np.pi/L  # Wavenumber "grid"
    xi = 2*np.pi*np.fft.fftfreq(nx+1, 2*L/(nx))
    
    # k = 2*np.pi*np.fft.rfftfreq(nx+1, 2*L/(nx))
    
    
    # (this is the order in which numpy's FFT gives the frequencies)

    # Initial data
    u = np.exp(-x**2)# np.sin(2*x)**2 * (x<-L/4)
    uhat0 = np.fft.fft(u)
    plt.plot(x,u)
    epsilon=b  # Diffusion coefficient
    # a = 1.0       # Advection coefficient

    # print(np.shape(xi))
    # print(np.shape(u))
    # Now we solve the problem
    d = np.ones(len(xi))*d
    dfft = np.fft.fft(d)
    # print(dfft)
    # print((dfft/((1.j*xi*a - epsilon*xi**2+c)+1.j*xi)))
    for i in range(0,N):
        t = i*dt
        al = uhat0 + dfft/((1.j*xi*a - epsilon*xi**2+c)+1.j*xi)
        uhat = np.exp((1.j*xi*a - epsilon*xi**2+c)*t) * al-(dfft/((1.j*xi*a - epsilon*xi**2+c)+1.j*xi))
        u = np.real(np.fft.ifft(uhat)) 
        
        U_psk[:, i] = u
        error =  U_CDI[:, i] - u
        error_2norm_CDI_PSTDK[i] =np.linalg.norm(error)/np.sqrt(nx)
        error =  U_ps[:, i] - u
        error_2norm_PSTDI_PSTDK[i] =np.linalg.norm(error)/np.sqrt(nx)

    
    plt.figure()  
    plt.imshow((U_psk), cmap='viridis',  extent=[0, T, 0, 2*L], aspect='auto')
    plt.colorbar()
    plt.title("CD Scheme with PSTDK")
    plt.xlabel('Time')
    plt.ylabel('Position')
    title = "PSTDK - nx = "+str(nx+1)
    plt.title(title)
    plt.grid(True)
    plt.show()
    plt.close()
    
    plt.figure();
    # plt.plot(x, ui, '-', color='g', label = 'Initial'); 
    plt.plot(x, U_CDI[:, int(nt/2)], '-', color='b', label = 'CDI'); 
    plt.plot(x, U_ps[:, int(nt/2)], '-', color = 'r', label = 'PSI');
    plt.plot(x, U_psk[:, int(nt/2)+1], '-', color = 'g', label = 'KPS');
    
    res =  [U_CDI[:, nt-1], U_ps[:, nt-1], U_psk[:, nt]]
    
    plt.legend(loc="upper right")
    title = "Comparison - nx = "+str(nx+1)
    plt.title(title)
    plt.show()
    plt.close()
    plt.figure()
    plt.plot(error_2norm_CDI_PSTD[0:nt], '.', label = 'PSTDI_CDI')
    
    plt.plot(error_2norm_CDI_PSTDK[0:nt], '.', label = 'PSTDK_CDI')
    
    plt.plot(error_2norm_PSTDI_PSTDK[0:nt], '.', label = 'PSTDI_PSTDK')
    plt.legend()
    plt.show()
    plt.close()
    return res
In [142]:
a = 1.0
b = 0.01
c = 0.1; 
d = 0.5

test_PS(a, b, c, d, 63, 1000)
test_PS(a, b, c, d, 127, 1000)
res = test_PS(a, b, c, d, 255, 1000)

plt.figure();
plt.plot(res[0], '-', color='b', label = 'CDI Method'); 
plt.plot(res[1], '-', color = 'r', label = 'Spectral Method');
plt.plot(res[2], '-', color = 'g', label = 'Ketchinson Spectral Method');
plt.legend(loc="upper right")
title = "Comparison - nx = "+str(255+1)
plt.title(title)
plt.show()
plt.close()
In [143]:
a = 1.0
b = 1.0
c = 0.1; 
d = 0.5

test_PS(a, b, c, d, 63, 1000)
test_PS(a, b, c, d, 127, 1000)
res = test_PS(a, b, c, d, 255, 1000)

plt.figure();
plt.plot(res[0], '-', color='b', label = 'CDI Method'); 
plt.plot(res[1], '-', color = 'r', label = 'Spectral Method');
plt.plot(res[2], '-', color = 'g', label = 'Ketchinson Spectral Method');
plt.legend(loc="upper right")
title = "Comparison - nx = "+str(255+1)
plt.title(title)
plt.show()
plt.close()
In [144]:
a = 0.1
b = 1.0
c = 0.1; 
d = 0.5

test_PS(a, b, c, d, 63, 1000)
test_PS(a, b, c, d, 127, 1000)
res = test_PS(a, b, c, d, 255, 1000)

plt.figure();
plt.plot(res[0], '-', color='b', label = 'CDI Method'); 
plt.plot(res[1], '-', color = 'r', label = 'Spectral Method');
plt.plot(res[2], '-', color = 'g', label = 'Ketchinson Spectral Method');
plt.legend(loc="upper right")
title = "Comparison - nx = "+str(255+1)
plt.title(title)
plt.show()
plt.close()

Trying to compare with the semi-infinite solutions of Ogata and Banks

Solution for No Production Decay

Lapidus and Amundson 1952
Ogata and Banks 1961
Equation: R*dC/dt = D*d2C/dx2 - v*dC/dx
IC: C(x, 0) = Ci
BCs:Dirichlet for 0<t<to: C(0, t) = Co
Neumann: dC/dx(inf, t) = 0

$$ \frac {\partial U}{\partial t} = D\frac {\partial^2 U}{\partial x^2} - v\frac {\partial U}{\partial x} + c*U + d $$$$ IC: U(x, 0) = C_i $$

$$ BC - Left: U(0, t) = C_o \ \ \ \ 0<t<t_o $$ $$ BC - Right: \frac {dU(\infty, t)}{dx} = 0 $$ $$ c = 0 $$ $$ d = 0 $$

In [26]:
def Analytical_Ogata(a, b, c, d, L, T, nx, nt):
    # Analytical - du/dt = -v*du/dx+D*d2u/dx2

    Ci = 0; 
    Co = 2;
    # dCinf/dx = 0;
    R = 1; 
    D = b
    v = -a
    print(v, D, c, d)
    x = np.linspace(0, L, nx+1);
    t = np.linspace(0, T, nt);


    C = np.zeros((nt, nx))
    for i in range(nt):  
        for j in range(nx):      
            C[i, j] = Ci + ((Co-Ci)/2)*(erfc((x[j]-v*t[i])/(2*np.sqrt(D*t[i]/R))) + (np.exp(v*x[j]/(D/R)))*(erfc((x[j]+v*t[i])/(2*np.sqrt(D*t[i]/R)))))

    U_exact = np.transpose(C);

    fig = plt.figure()  
    plt.imshow(np.transpose(C), cmap='viridis',  extent=[0, T, L, 0], aspect='auto')
    plt.colorbar()
    plt.title("Exact")
    plt.xlabel('Time')
    plt.ylabel('Position')

    #fig.set_size_inches(30.,18.)
    # plt.savefig('U_CDI.png', dpi = 900)
    plt.show()
    
    return U_exact[:, -1]
In [27]:
## Ogata and Banks - No production or decay terms
axs = [-1.0, -0.1, -0.01]
bs = [0.01, 0.1]
it = 0
L = 1.0;
T = 5.0;
nx = 128;
nt = 5001; 
plt.figure()
print("ddfed")
U_compab = np.zeros((12, nx))
for a in axs:
    for b in bs: 
        pe = abs(a/b);
        # U_exact =  
        U_compab[it, :]  = Analytical_Ogata(a, b, 0, 0, L, T, nx, nt)
        U_compab[it+1, :] = testCD_PS(-a, b, 0, 0, nx-1, nt-1)
        # U_comp = comparePeclet(a, b)
        
        plt.plot(U_compab[i, :],'+-'); plt.plot(U_compab[i+1, :],'-', label = str(np.round(pe, 3)));
        plt.ylim([0, 2.1])
        # plt.title("Pe = "+str(np.round(pe, 3)))
        
        plt.xlabel("Dimensionless Time")
        plt.ylabel("Dimensionless Concentration")
        plt.title("Ogata and Banks - Different Pe="+str(np.round(pe, 3)))
        
        it = it+2

# Printing the temporal grid tests
zipped1 = list(zip(U_compab[0, :], U_compab[1, :], U_compab[2, :], U_compab[3, :], U_compab[4, :], U_compab[5, :], U_compab[6, :], U_compab[7, :], U_compab[8, :], U_compab[9, :], U_compab[10, :], U_compab[11, :]))

df2 = pd.DataFrame(zipped1, columns=['OB', 'Pe = 100.0', 'OB', 'Pe = 10.0', 'OB', 'Pe = 10.0','OB', 'Pe = 1.0','OB', 'Pe = 1.0','OB', 'Pe = 0.1'])

df2.to_csv('Analysis_OgataBanks_complete.csv', index=False)
ddfed
1.0 0.01 0 0
C:\Users\gaura\anaconda3\envs\newenv\lib\site-packages\ipykernel_launcher.py:18: RuntimeWarning: invalid value encountered in double_scalars
C:\Users\gaura\anaconda3\envs\newenv\lib\site-packages\ipykernel_launcher.py:18: RuntimeWarning: divide by zero encountered in double_scalars
<Figure size 432x288 with 0 Axes>
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-27-e5ffb99cfc65> in <module>
     15         # U_exact =
     16         U_compab[it, :]  = Analytical_Ogata(a, b, 0, 0, L, T, nx, nt)
---> 17         U_compab[it+1, :] = testCD_PS(-a, b, 0, 0, nx-1, nt-1)
     18         # U_comp = comparePeclet(a, b)
     19 

ValueError: cannot copy sequence with size 2 to array axis with dimension 128

Didn't work - Using MMS

MMS for the same manufactured solution as it satisfies the present condition

For an equation with given IC and periodic BCs : $$ \frac {\partial U}{\partial t} = D\frac {\partial^2 U}{\partial x^2} + v\frac {\partial U}{\partial x}$$ $$ IC: U(x, 0) = C_i $$ $$ BC - Left: U(0, t) = -C_o $$ $$ BC - Right: U(0, t) = -C_o $$ $$ BC - Left: \frac {dU(0, t)}{dx} = 0 $$ $$ BC - Right: \frac {dU(1, t)}{dx} = 0 $$

The chosen equation is: $$ U(x,t) = C_o \cos(\pi(2x -1)) e^{-t}$$

The initial and boundary condition terms: I.C.: Cisin(pix) Boundary Conditions: Periodic BC $$ U(x, 0) = C_o\cos( \pi(2x -1))$$ $$ U(0, t) = -C_o e^{-t} $$ $$ U(1, t) = -C_o e^{-t} $$

$$ \frac {dU}{dx} (0, t) = -2\pi C_o\sin( \pi(2*0 -1)) e^{-t} = 0$$

$$ \frac {dU}{dx} (1, t) = -2\pi C_o\sin(\pi(2*1 -1)) e^{-t} = 0$$

The derivative terms are:

$$ \frac {dU}{dt} = -C_o\cos(\pi (2*x-1)) e^{-t}$$$$ \frac {dU}{dx} = -C_o 2 \pi \sin(\pi (2*x-1)) e^{-t}$$$$ \frac {d^2U}{dx^2} = -C_o 4 \pi^2 \cos(\pi (2*x-1)) e^{-t}$$

The new source term will be:

$$ d = \frac {dU}{dt} - v*\frac {dU}{dx} - D*\frac {d^2U}{dx^2}$$$$ d = -C_o\cos(\pi (2*x-1)) e^{-t} - v*2*C_o\pi \sin(\pi (2*x-1)) e^{-t} - D*C_o*4*\pi^2 \cos(\pi (2*x-1)) e^{-t}$$$$ d = -C_o e^{-t}(\cos(\pi (2*x-1)) - v*2*\pi \sin(\pi (2*x-1)) - D*4*\pi^2 \cos(\pi (2*x-1))) $$
In [28]:
def numerical_scheme_CDI(u_CDI, a, b, c, d, dx, dt):
        alpha = b*dt/dx
        beta = a*dt/dx**2
        nx = len(u_CDI)-1;
        A = np.zeros((nx+1, nx+1))
        bi = d*dt+u_CDI;
        for j in range(1, len(u_CDI)-1):

            A[j, j-1] = -(beta-alpha/2)
            A[j, j] = 1+2*beta-c*dt
            A[j, j+1] = -(beta+alpha/2)

        # Applying first order Periodic Boundary Condition
        A[0, 0]=(1+2*beta)-c*dt;
        A[0, 1]= (-alpha/2-beta)
        A[0, nx]= (alpha/2-beta)    


        A[nx, 0]= (-alpha/2-beta) ;
        A[nx, nx-1]= (alpha/2 - beta)
        A[nx, nx]= 1+2*beta-c*dt    


        
        u_CDI = solve(A, bi)

        return u_CDI
In [29]:
def numerical_scheme_PSTD(u, a, b, c, d, k, dx, dt):
        uk = np.fft.rfft(u)
        
        uk[:] = (uk[:])/(1-(-a*k**2 + b*1j*k+c)*dt) 
        u = np.fft.irfft(uk) + d*dt
        return u
In [30]:
# Creating function for running the code and getting the error with different number of nodes with left dirichlet and right neumann condition. 
def numericalcodes(nx, nt): 
    # Define the grid
    a = 1.0;
    b = 0.01;
    c = 0.0;
    Co = 1.0;
    L = 1.0
    T = 1.0
    dm = np.zeros((nx+1, nt+1));
    U_exact = np.zeros((nx+1, nt+1));
    # Define x & t
    x = np.linspace(0, L, nx+1)
    t = np.linspace(0, T, nt+1)
    for i in range(nt+1):
        for j in range (nx+1):
            dm[j, i] = -(Co*np.exp(-t[i]))*(np.cos(np.pi*(2*x[j]-1)) - 2*b*np.pi*np.sin(np.pi*(2*x[j]-1)) - a*4*(np.pi**2)*np.cos(np.pi*(2*x[j]-1)))
    
    for i in range(nt+1):
        for j in range (nx+1):
            U_exact[j, i] = Co*np.cos(np.pi*(2*x[j]-1))*np.exp(-t[i])
    fig = plt.figure()  
    plt.imshow((U_exact), cmap='viridis',  extent=[0, L, 0, T], aspect='auto')
    plt.colorbar()
    plt.title("Exact")
    plt.xlabel('Time')
    plt.ylabel('Position')
    plt.show()

    # Defining variables and errors
    U_CDI = np.zeros((nx+1, nt+1))
    U_ps = np.zeros((nx+1, nt+1))
    error_2norm_CDI_PSTD = np.zeros(nt+1)
    error_2norm_Exact_CDI = np.zeros(nt+1)
    error_2norm_Exact_PSTD = np.zeros(nt+1)
    
    # Discretization of space and time
    dx = L/(nx)
    dt = T/(nt)
    
    print("\n here \n")
    print(dx)
    print(x[1])
    print(dt)
    print(t[1])
    print("\n here \n")
    u = Co*np.cos(np.pi*(2*x-1)) # np.zeros(nx);
    u_ps = u;
    u_CDI=u;
    u1 = u
    k = 2*np.pi*np.fft.rfftfreq(nx+1, L/(nx))    
    
    # Solving over time
    for i in range(nt+1):
        u_CDI = numerical_scheme_CDI(u_CDI, a, b, c, dm[:, i], dx, dt)
        u_ps = numerical_scheme_PSTD(u_ps, a, b, c, dm[:, i], k, dx, dt)
        if(i<2):
            u1= u_ps
        U_CDI[:, i] = u_CDI
        U_ps[:, i] = u_ps
        error =  U_CDI[:, i] - U_ps[:, i]
        error_2norm_CDI_PSTD[i] =np.linalg.norm(error)/np.sqrt(nx)
        
        error =  U_ps[:, i] - U_exact[:, i]
        error_2norm_Exact_PSTD[i] =np.linalg.norm(error)/np.sqrt(nx)
        errorPSTD = error_2norm_Exact_PSTD[i]
        
        error =  U_CDI[:, i] - U_exact[:, i]
        error_2norm_Exact_CDI[i] =np.linalg.norm(error)/np.sqrt(nx)
        errorCDI = error_2norm_Exact_CDI[i]

    # Errors at time 1s for the two methods(CDI, PSTD) and 7 number of nodes (dx in 20, 40, 60, 80, 100, 120, 150)
    errornorms = np.ones(7)
    errornorms[0] = errorCDI
    errornorms[1] = errorPSTD
    
    plt.rcParams['font.size'] = 18
    fig = plt.figure() 
    plt.imshow(U_CDI, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
    plt.colorbar()
    title = "CDI - nx = "+str(nx+1)
    plt.title(title)
    plt.xlabel('Time')
    plt.ylabel('Position')
    plt.show()

    plt.figure()  
    plt.imshow((U_ps), cmap='viridis',  extent=[0, T, 0, L], aspect='auto')
    plt.colorbar()
    plt.title("CD Scheme with PSTD")
    plt.xlabel('Time')
    plt.ylabel('Position')
    title = "PSTDI - nx = "+str(nx+1)
    plt.title(title)
    plt.grid(True)
    plt.show()
    plt.close()

    plt.figure();
    plt.plot(x, u, '-', color='g', label = 'Initial'); 
    plt.plot(x, u1, '-', color='c', label = 'InitialPS'); 
    
    plt.plot(x, U_CDI[:, nt-1], '-*', color='b', label = 'CDI Method'); 
    plt.plot(x, U_ps[:, nt-1], '-+', color = 'r', label = 'Spectral Method');
    plt.plot(x, U_exact[:, nt-1], '.-', color = 'k', label = 'Exact');
    plt.legend(loc="upper right")
    title = "Comparison - nx = "+str(nx+1)
    plt.title(title)
    plt.show()
    plt.close()
    
    plt.figure()
    plt.plot(error_2norm_CDI_PSTD[0:nt+1], '.', label = 'PSTDI_CDI')
    plt.plot(error_2norm_Exact_PSTD[0:nt+1], '.', label = 'Exact_PSTD')
    plt.plot(error_2norm_Exact_CDI[0:nt+1], '.', label = 'Exact_CDI')
    plt.legend()
    plt.show()
    plt.close()
    return errornorms
    
In [31]:
# Equation with decay terms 

# Spatial MMS tests
errorsx = np.ones((8, 2))
errornorms = np.ones(8)
nx = 32;
nt = 1000;
nxs = [15, 31, 63, 127, 255, 511, 1023, 2047]
nts = [250, 500, 1000, 2000, 4000]
for i in range(6):
    errornorms = numericalcodes(nxs[i], nt)
    errorsx[i, 0] = errornorms[0]
    errorsx[i, 1] = errornorms[1]
fig = plt.figure()
plt.figure()
fig.set_size_inches(10.,2.)
plt.plot(nxs, errorsx[:, 0], '-', label = 'CDI')
plt.plot(nxs, errorsx[:, 1], '.-', label = 'PSTD')

plt.xlabel('Grid Numbers')
plt.ylabel('Error Norms')
plt.legend()

fig = plt.figure()
fig.set_size_inches(10.,7.)
px = np.zeros((6, 2))
for g in range(6):
    for k in range(2):
        px[g, k] = np.log((np.abs(errorsx[g+2, k]) - np.abs(errorsx[g+1, k]))/(np.abs(errorsx[g+1, k]) - np.abs(errorsx[g, k])))/(np.log(nxs[g]/nxs[g+1]))

plt.plot(nxs[0:6], px[:, 0], '-', label = 'CDI')
plt.plot(nxs[0:6], px[:, 1], '-.', label = 'PSTD')

plt.xlabel('Spatial Grid Numbers')
plt.ylabel('Observed Order')
plt.legend()

# Printing the spatial grid tests

zipped1 = list(zip(nxs, errorsx[:, 0], errorsx[:, 1], px[:, 0], px[:, 1]))
df2 = pd.DataFrame(zipped1, columns=['Spatial Discretization Numbers', 'CDI_Error', 'PSTD_Error', 'CDI_Order', 'PSTD_Order'])

df2.to_csv('Analysis_Spatial25_Complete.csv', index=False)
 here 

0.06666666666666667
0.06666666666666667
0.001
0.001

 here 

 here 

0.03225806451612903
0.03225806451612903
0.001
0.001

 here 

 here 

0.015873015873015872
0.015873015873015872
0.001
0.001

 here 

 here 

0.007874015748031496
0.007874015748031496
0.001
0.001

 here 

 here 

0.00392156862745098
0.00392156862745098
0.001
0.001

 here 

 here 

0.0019569471624266144
0.0019569471624266144
0.001
0.001

 here 

C:\Users\gaura\anaconda3\envs\newenv\lib\site-packages\ipykernel_launcher.py:29: RuntimeWarning: invalid value encountered in log
C:\Users\gaura\anaconda3\envs\newenv\lib\site-packages\ipykernel_launcher.py:29: RuntimeWarning: divide by zero encountered in log
<Figure size 720x144 with 0 Axes>
In [32]:
a = 1.0;
b = 0.01;
c = 0.0;
Co = 1.0;
L = 1.0
T = 1.0
dm = np.zeros((nx+1, nt+1));
U_exact = np.zeros((nx+1, nt+1));
# Define x & t
x = np.linspace(0, L, nx+1)
t = np.linspace(0, T, nt+1)
for i in range(nt+1):
    for j in range (nx+1):
        U_exact[j, i] = Co*np.cos(np.pi*(2*x[j]-1))*np.exp(-t[i])


plt.plot(x, U_exact[:, 1], '.-', color = 'k', label = 'Exact0');
plt.plot(x, U_exact[:, 500], '.-', color = 'g', label = 'Exact500');
plt.plot(x, U_exact[:, 900], '.-', color = 'b', label = 'Exact900');
In [33]:
def tests2(a, b, c, d, nx, nt):

    L = 1.0
    T = 1.0


    x = np.linspace(0, L, nx+1)
    t = np.linspace(0, T, nt+1)
    xin = np.linspace(0, L, nx+1);
    u = np.exp(-x**2)#np.sin(np.pi*xin)
    ui = u;
    U_CDI = np.zeros((nx+1, nt+1))
    error_2norm_CDI_PSTD = np.zeros(nt+1)

    # Discretization of space and time
    dx = L/(nx-1)
    dt = T/(nt+1)


    #plt.figure()
    #plt.plot(u, label='Initial')
    def numerical_scheme_CDI(u_CDI, a, b, c, d, dx, dt):
        alpha = a*dt/dx
        beta = b*dt/dx**2
        nx = len(u_CDI)-1;
        A = np.zeros((nx+1, nx+1))
        bi = d*dt+u_CDI;
        for j in range(1, len(u_CDI)-1):

            A[j, j-1] = -(beta[j]-alpha[j]/2)
            A[j, j] = 1+2*beta[j]-c[j]*dt
            A[j, j+1] = -(beta[j]+alpha[j]/2)

        # Applying first order Periodic Boundary Condition
        A[0, 0]=(1+2*beta[0])-c[0]*dt;
        A[0, 1]= (-alpha[0]/2-beta[0])
        A[0, nx]= (alpha[0]/2-beta[0])    


        A[nx, 0]= (-alpha[nx]/2-beta[nx]) ;
        A[nx, nx-1]= (alpha[nx]/2 - beta[nx])
        A[nx, nx]= 1+2*beta[nx]-c[nx]*dt    

        u_CDI = solve(A, bi)

        return u_CDI



    # Apply the numerical scheme
    U_CDI = np.zeros((nx+1, nt+1))
    error_2norm_QI_PSTD = np.zeros(nt)

    u_CDI=u;

    for n in range(nt):
        u_CDI = numerical_scheme_CDI(u_CDI, a, b, c, d, dx, dt)
        U_CDI[:, n] = u_CDI

    plt.rcParams['font.size'] = 18
    fig = plt.figure() 
    plt.imshow(U_CDI, cmap='viridis', extent=[0, T, 0, L], aspect='auto')
    plt.colorbar()
    title = "CDI - nx = "+str(nx+1)
    plt.title(title)
    plt.xlabel('Time')
    plt.ylabel('Position')
    plt.show()


    k = 2*np.pi*np.fft.rfftfreq(nx+1, L/(nx))
    k2=k**2;

    # Defining variables
    U_ps = np.zeros((nx+1, nt+1))

    # Solving over time
    for i in range(nt+1): # 
        
        uk = np.fft.rfft(u)

        uk[:] = (uk[:])/(1-(-b*k**2 + a*1j*k+c)*dt) 
        u = np.fft.irfft(uk) + d*dt
        U_ps[:, i] = u
        error =  U_CDI[:, i] - u
        error_2norm_CDI_PSTD[i] =np.linalg.norm(error)/np.sqrt(nx)
 

    fig = plt.figure()  
    plt.imshow((U_ps), cmap='viridis',  extent=[0, T, 0, L], aspect='auto')
    plt.colorbar()
    plt.title("CD Scheme with PSTD")
    plt.xlabel('Time')
    plt.ylabel('Position')
    title = "PSTDI - nx = "+str(nx+1)
    plt.title(title)
    plt.grid(True)
    plt.show()
    plt.close()

    plt.figure();
    plt.plot(x, ui, '-', color='g', label = 'Initial'); 

    plt.plot(x, U_CDI[:, nt-1], '-*', color='b', label = 'CDI Method'); 
    plt.plot(x, U_ps[:, nt-1], '-+', color = 'r', label = 'Spectral Method');
    res =  [U_CDI[:, nt-1], U_ps[:, nt-1]]
    # print(np.shape(res))
    plt.legend(loc="upper right")
    title = "Comparison - nx = "+str(nx+1)
    plt.title(title)
    plt.show()
    plt.close()
    plt.figure()
    plt.plot(error_2norm_CDI_PSTD[0:nt], '.', label = 'PSTDI_CDI')
    plt.legend()
    plt.show()
    plt.close()
    return res
In [34]:
nx = 63
x = np.linspace(0,1, nx+1)
a = 1.0*np.sin(np.pi*x)
b = 0.01*np.cos(np.pi*x)
c = a+10*b; 
d = a-10*b;

tests2(a, b, c, d, nx, 1000)

nx = 127
x = np.linspace(0,1, nx+1)
a = 1.0*np.sin(np.pi*x)
b = 0.01*np.cos(np.pi*x)
c = a+10*b; 
d = a-10*b;
tests2(a, b, c, d, nx, 1000)

nx = 255
x = np.linspace(0,1, nx+1)
a = 1.0*np.sin(np.pi*x)
b = 0.01*np.cos(np.pi*x)
c = a+10*b; 
d = a-10*b;
res = tests2(a, b, c, d, nx, 1000)
print('ldd')
plt.figure();
plt.plot(res[0], '-*', color='b', label = 'CDI Method'); 
plt.plot(res[1], '-+', color = 'r', label = 'Spectral Method');
plt.legend(loc="upper right")
title = "Comparison - nx = "+str(nx+1)
plt.title(title)
plt.show()
plt.close()
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-34-dcc45f02d821> in <module>
      6 d = a-10*b;
      7 
----> 8 tests2(a, b, c, d, nx, 1000)
      9 
     10 nx = 127

<ipython-input-33-b7d3dfc3a514> in tests2(a, b, c, d, nx, nt)
     80         uk = np.fft.rfft(u)
     81 
---> 82         uk[:] = (uk[:])/(1-(-b*k**2 + a*1j*k+c)*dt)
     83         u = np.fft.irfft(uk) + d*dt
     84         U_ps[:, i] = u

ValueError: operands could not be broadcast together with shapes (64,) (33,) 
In [ ]:
 
In [ ]:
 
In [ ]:
 
In [ ]: